End of lab session

frequencySpectrum.m

@@ -0,0 +1,60 @@
+function [power, duration] = frequencySpectrum(signal, fs, pad)
+%%%%%%%%%%%%%%%%%%
+%function power = frequencySpectrum(signal, fs, pad)
+%
+% Task: Display the power spectrum (lin and log scale) of a given signal
+%
+% Input:
+%	- signal: the input signal to process
+%	- fs: the sampling rate
+%	-pad: boolean if true, signal is padded with 0 to the next power of 2 -> FFT instead of DFT
+%
+% Output: 
+%	- power: the power spectrum
+%	
+%
+% Guillaume Gibert, guillaume.gibert@ecam.fr
+% 25/04/2022
+%%%%%%%%%%%%%%%%%%
+
+n = length(signal);        % number of samples
+
+if (pad)
+	n = 2^nextpow2(n);
+end
+
+tic
+y = fft(signal, n);% compute DFT of input signal
+duration = toc;
+
+power = abs(y).^2/n;    % power of the DFT
+
+[val, ind] = max(power); % find the mx value of DFT and its index
+
+% plots
+figure;
+
+subplot(1,3,1) % time plot 
+t=0:1/fs:(n-1)/fs; % time range
+%pad signal with zeros
+if (pad)
+	signal = [ signal; zeros( n-length(signal), 1)];
+end
+plot(t, signal)
+xticks(0:0.1*fs:n*fs);
+xticklabels(0:0.1:n/fs);
+xlabel('Time (s)');
+ylabel('Amplitude (a.u.)');
+
+subplot(1,3,2) % linear frequency plot
+f = (0:n-1)*(fs/n);     % frequency range
+plot(f,power, 'b*'); hold on;
+plot(f,power, 'r');
+xlabel('Frequency (Hz)')
+ylabel('Power (a.u.)')
+
+subplot(1,3,3) % log frequency plot
+plot(f,10*log10(power/power(ind)));
+xlabel('Frequency (Hz)')
+ylabel('Power (dB)')
+

spectrogram.m

@@ -0,0 +1,38 @@
+function spectrogram(signal, samplingFreq, step_size, window_size)
+%%%%%%%%%%%%%%%%%%%%%%%
+%function spectrogram(signal, samplingFreq, step_size, window_size)
+% ex.:  spectrogram(signal, samplingFreq, step_size, window_size)
+%
+% Task: Plot the spectrogram of a given signal
+%
+% Inputs:
+%	-signal: temporal signal to analyse 
+%	-samplingFreq: sampling frequency of the temporal signal
+% 	-step_size: how often the power spectrum will be computed in ms
+%	-window_size: size of the analysing window in ms
+%
+% Ouput: None
+%
+% author: Guillaume Gibert (guillaume.gibert@ecam.fr)
+% date: 14/03/2023
+%%%%%%%%%%%%%%%%%%%%%%%
+
+figure;
+	subplot(2,1,1);
+t=0:1/samplingFreq:length(signal)/samplingFreq-1/samplingFreq;
+plot(t, signal');
+xlim([0 length(signal)/samplingFreq-1/samplingFreq]);
+ylabel('amplitude (norm. unit)');
+	subplot(2,1,2);
+step = fix(step_size*samplingFreq/1000);     % one spectral slice every step_size ms
+window = fix(window_size*samplingFreq/1000);  % window_size ms data window
+fftn = 2^nextpow2(window); % next highest power of 2
+[S, f, t] = specgram(signal, fftn, samplingFreq, window, window-step);
+S = abs(S(2:fftn*4000/samplingFreq,:)); % magnitude in range 0<f<=4000 Hz.
+S = S/max(S(:));           % normalize magnitude so that max is 0 dB.
+S = max(S, 10^(-40/10));   % clip below -40 dB.
+S = min(S, 10^(-3/10));    % clip above -3 dB.
+imagesc (t, f, log(S));    % display in log scale
+set (gca, "ydir", "normal"); % put the 'y' direction in the correct direction
+xlabel('time (s)');
+ylabel('frequency (Hz)');

speech_analysis.m

@@ -1 +1,131 @@
-pkg load signal
+% 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);