% Load the speech signal
[signal, fs] = audioread('modulator22.wav');
t=0:1/fs:length(signal)/fs-1/fs;
figure;
plot(t, signal);
xlabel('Time(s)');
ylabel('Amplitude(n.u.)');
% Modify the sampling frequency
% new_fs = fs/2;
% 
% % Save the modified signal
% audiowrite('modified_modulator22.wav', signal, new_fs);
% 
% % Listen to the generated sound
% [signal_modified, fs_modified] = audioread('modified_modulator22.wav');
% sound(signal_modified, fs_modified);
% figure;
% plot(t, signal_modified);
% xlabel('Time(s)');
% ylabel('Amplitude(n.u.)');
% title('Signal modified (fs/2)');

% frequencySpectrum(signal, fs, 0);
% frequencySpectrum(signal, fs, 1);

% Number of repetitions for measurement
% num_repetitions = 10;
% fft_times = zeros(num_repetitions, 1);
% dft_times = zeros(num_repetitions, 1);
% 
% for i = 1:num_repetitions
%     % Measure time to compute FFT
%     tic;
%     Y = frequencySpectrum(signal, fs, 0);
%     fft_times(i) = toc;
% 
%     % Measure time to compute DFT
%     tic;
%     Y_dft = frequencySpectrum(signal, fs, 1);
%     dft_times(i) = toc;
% end
% 
% % Display results
% disp('FFT computation times:');
% disp(fft_times);
% disp(['Average FFT time: ', num2str(mean(fft_times)), ' seconds']);
% disp(['Standard deviation of FFT time: ', num2str(std(fft_times)), ' seconds']);
% 
% disp('DFT computation times:');
% disp(dft_times);
% disp(['Average DFT time: ', num2str(mean(dft_times)), ' seconds']);
% disp(['Standard deviation of DFT time: ', num2str(std(dft_times)), ' seconds']);

% spectrogram(signal, fs, 5, 30);
% spectrogram(signal, fs, 5, 5);


t1 = 0.8577;
t2 = 0.9720;
t3 = 1.4090;
t4 = 1.6734;
t5 = 1.9871;
t6 = 2.2282;

% Extract the portion of the signal corresponding to the time interval [t1, t2]
start_sample1 = round(t1 * fs);
end_sample1 = round(t2 * fs);
start_sample2 = round(t3 * fs);
end_sample2 = round(t4 * fs);
start_sample3 = round(t5 * fs);
end_sample3 = round(t6 * fs);
signal_vowel1 = signal(start_sample1:end_sample1, 1);
signal_vowel2 = signal(start_sample2:end_sample2, 1);
signal_vowel3 = signal(start_sample3:end_sample3, 1);

% Generate the time vector corresponding to the extracted portion
t_vowel1 = (start_sample1:end_sample1) / fs;
t_vowel2 = (start_sample2:end_sample2) / fs;
t_vowel3 = (start_sample3:end_sample3) / fs;

% Plot the extracted portion of the signal
figure;
subplot(1,3,1)
plot(t_vowel1, signal_vowel1);
xlabel('Time (s)');
ylabel('Amplitude (n.u.)');
title('Portion Signal /ʌ/');

subplot(1,3,2)
plot(t_vowel2, signal_vowel2);
xlabel('Time (s)');
ylabel('Amplitude (n.u.)');
title('Portion Signal /uː/');

subplot(1,3,3)
plot(t_vowel3, signal_vowel3);
xlabel('Time (s)');
ylabel('Amplitude (n.u.)');
title('Portion Signal /iː/');

% Compute the power spectrum of each vowel signal
P_vowel1 = frequencySpectrum(signal_vowel1, fs, 0);
P_vowel2 = frequencySpectrum(signal_vowel2, fs, 0);
P_vowel3 = frequencySpectrum(signal_vowel3, fs, 0);

% Apply a low-pass filter to retrieve the envelope
N=8;
fc = 1000; % Cutoff frequency for the low-pass filter
[b, a] = butter(N, fc / (fs / 2));
% freqz(b,a);
% Z=roots(b);
% P=roots(a);
% figure;
% zplane(Z, P);
% title('Zeros and poles of the transfer function of the IIR filter');
% legend('zeros', 'poles');
% grid on
% filter the signal
signal_filtered=filter(b, a, signal);
figure;
plot(signal_filtered);

envelope1 = filter(b, a, abs(P_vowel1));
figure;
plot(envelope1);
envelope2 = filter(b, a, abs(P_vowel2));
figure;
plot(envelope2);
envelope3 = filter(b, a, abs(P_vowel3));
figure;
plot(envelope3);