diff --git a/iirFilter.m b/iirFilter.m deleted file mode 100644 index 36701b2..0000000 --- a/iirFilter.m +++ /dev/null @@ -1,37 +0,0 @@ -function [Z, P]= iirFilter(N, cutoffFreq, samplingFreq, filterType) -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -% function [Z, P] = iirFilter(N, cutoffFreq, samplingFreq, filterType) -% ex.: [Z, P] =iirFilter(6, 10, 500, 1) -% -% Task: To create and analyze an IIR low pass filter (Butterworth ror Chebychev) -% -% Inputs: -% -N: order of the filter -% -cutoffFreq: below this frequency, signal is not modified and above signal is attenuated -% -samplingFreq: sampling frequency (In Hz) -% -filterType: Butterworth if equal to 1 and Chebychev if equal to 2 -% -% Outputs: -% -% -% Author: Guillaume Gibert, guillaume.gibert@ecam.fr -% Date: 09/04/2025 -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - -if (filterType == 1) - [b, a] = butter(N, cutoffFreq/(samplingFreq/2)); -elseif (filterType == 2) - Rp = 10; % bandpass ripple of Rp dB - [b, a] = cheby1(N, Rp, cutoffFreq/(samplingFreq/2)); -else - disp('Filter type is incorrect!') - return -end - -[Z, P] = zeroPole(a, b, 1); - -figure; -freqz(b, a, N, samplingFreq); -title('Frequency response'); -grid on; - diff --git a/speech_analysis.m b/speech_analysis.m index 96651e0..f2cbe87 100644 --- a/speech_analysis.m +++ b/speech_analysis.m @@ -35,7 +35,7 @@ title(['Temporal Variation of ', filepath]); grid on; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -%Frequency Spectrum +%% Frequency Spectrum %FFT tic; [yFFT, FFT_Time]=frequencySpectrum(y,Fs, 1); @@ -47,7 +47,7 @@ disp(DFT_Time); %Modify the padding to make the change. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - +%% Spectrogram spectrogram(y, Fs, 5,50) title('Spectrogram of modulator22.wav'); colorbar; @@ -61,6 +61,7 @@ xlabel('Time (s)'); %F2 : 1695.8136433413672, 1550.9109531347972, 566.7831612330604, %1721.8044733141373, 1802.7920754749957, 1891.9059418088873 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%% First downsampling (Shannon-Nyquist problem) desiredFreq = 4000; %in Hz % --- Downsampling using downsample() --- @@ -79,7 +80,14 @@ disp(['--- Downsampling using decimate() ---']); disp(['New sampling frequency (decimate): ', num2str(Fs_decimated), ' Hz']); disp(['Number of samples (decimate): ', num2str(length(y_decimated))]); -% --- Plotting Downsampled Signals --- +%% --- Plotting Downsampled Signals --- +figure; +subplot(3,1,1); +plot(t, y); +xlabel('Time (seconds)'); +ylabel('Amplitude'); +title(['Original Signal (Fs = ', num2str(Fs), ' Hz)']); +grid on; t_ds = (0:length(y_downsampled_ds)-1) / Fs_downsampled_ds; subplot(3,1,2); plot(t_ds, y_downsampled_ds); @@ -87,7 +95,6 @@ xlabel('Time (seconds)'); ylabel('Amplitude'); title(['Downsampled Signal (downsample, Fs = ', num2str(Fs_downsampled_ds), ' Hz)']); grid on; - t_dec = (0:length(y_decimated)-1) / Fs_decimated; subplot(3,1,3); plot(t_dec, y_decimated); @@ -95,20 +102,18 @@ xlabel('Time (seconds)'); ylabel('Amplitude'); title(['Decimated Signal (decimate, Fs = ', num2str(Fs_decimated), ' Hz)']); grid on; - %{ -% --- Frequency Spectrum of Downsampled Signals --- +%% --- Frequency Spectrum of Downsampled Signals --- figure; subplot(2,1,1); [yFFT_ds, FFT_Time_ds]=frequencySpectrum(y_downsampled_ds,Fs_downsampled_ds, 1); disp(['FFT Time (downsampled): ', num2str(FFT_Time_ds)]); plot(yFFT_ds, Fs_downsampled_ds); title('FFT of Downsampled Signal (downsample)'); - subplot(2,1,2); [yFFT_dec, FFT_Time_dec]=frequencySpectrum(y_decimated,Fs_decimated, 1); disp(['FFT Time (decimated): ', num2str(FFT_Time_dec)]); -plot(yFFT_dec, Fs_decimated) +plot(yFFT_dec, Fs_decimated) title('FFT of Decimated Signal (decimate)'); %} %{ @@ -120,7 +125,6 @@ title(['Spectrogram of Downsampled Signal (downsample, Fs = ', num2str(Fs_downsa colorbar; ylabel('Frequency (Hz)'); xlabel('Time (s)'); - subplot(2,1,2); spectrogram(y_decimated, round(0.02*Fs_decimated), round(0.01*Fs_decimated), 512, Fs_decimated, 'yaxis'); title(['Spectrogram of Decimated Signal (decimate, Fs = ', num2str(Fs_decimated), ' Hz)']); @@ -128,15 +132,86 @@ colorbar; ylabel('Frequency (Hz)'); xlabel('Time (s)'); %} +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +desiredFilterFreq = 1000; +%% --- Low-pass FIR filter --- +order_fir = 30; +normalized_cutoff_fir = desiredFilterFreq / (Fs / 2); +y_fir_coeffs = fir1(order_fir, normalized_cutoff_fir, 'low'); % 'low' specifies a low-pass filter +y_fir_filtered = filter(y_fir_coeffs, 1, y); % Apply the FIR filter +figure; +freqz(y_fir_coeffs, 1, 512, Fs); % Plot the frequency response of the FIR filter +title('Frequency Response of FIR Low-Pass Filter'); +% FIR filter stability check (always stable) +disp('--- FIR Filter Stability ---'); +disp('FIR filters designed using fir1 are inherently stable.'); +% --- Low-pass IIR filter (Butterworth) --- +order_iir = 8; +normalized_cutoff_iir = desiredFilterFreq / (Fs / 2); +[b_iir, a_iir] = butter(order_iir, normalized_cutoff_iir, 'low'); % 'low' specifies a low-pass filter +y_iir_filtered = filter(b_iir, a_iir, y); % Apply the IIR filter +figure; +freqz(b_iir, a_iir, 512, Fs); % Plot the frequency response of the IIR filter +title('Frequency Response of IIR (Butterworth) Low-Pass Filter'); +%{ +% IIR filter stability check +disp('--- IIR Filter (Butterworth) Stability ---'); +poles_iir = roots(a_iir); +magnitudes_iir = abs(poles_iir); +if all(magnitudes_iir < 1) + disp('The IIR (Butterworth) filter is stable (all poles are inside the unit circle).'); +else + disp('The IIR (Butterworth) filter is NOT stable (some poles are outside or on the unit circle).'); + disp('Poles magnitudes:'); + disp(magnitudes_iir); +end +%} + +%% --- Downsampling after filtering --- +downsample_factor_filtered = round(Fs / desiredFreq); +Fs_ds_filtered = Fs / downsample_factor_filtered; +%{ +y_ds_fir_filtered = downsample(y_fir_filtered, downsample_factor_filtered); +disp(['--- Downsampling FIR filtered signal using downsample() ---']); +disp(['New sampling frequency (FIR filtered, downsample): ', num2str(Fs_ds_filtered), ' Hz']); +disp(['Number of samples (FIR filtered, downsample): ', num2str(length(y_ds_fir_filtered))]); +%} +y_ds_iir_filtered = downsample(y_iir_filtered, downsample_factor_filtered); +disp(['--- Downsampling IIR filtered signal using downsample() ---']); +disp(['New sampling frequency (IIR filtered, downsample): ', num2str(Fs_ds_filtered), ' Hz']); +disp(['Number of samples (IIR filtered, downsample): ', num2str(length(y_ds_iir_filtered))]); + +%% Plotting the signals +% --- Comparing Output Signals --- +% Temporal Variation +figure; +subplot(3,1,1); +plot(t, y); +xlabel('Time (seconds)'); +ylabel('Amplitude'); +title('Original Signal'); +grid on; +subplot(3,1,2); +plot(t, y_fir_filtered); +xlabel('Time (seconds)'); +ylabel('Amplitude'); +title('FIR Filtered Signal'); +grid on; +subplot(3,1,3); +plot(t, y_iir_filtered); +xlabel('Time (seconds)'); +ylabel('Amplitude'); +title('IIR Filtered Signal'); +grid on; % Play audios (using the audio data 'y' and its sampling rate 'Fs') %sound(y, Fs); % Play the original sound %sound(y, Fs*2); %sound(y_decimated,Fs_decimated) -sound(y_downsampled_ds,Fs_downsampled_ds) %Has distortion. This is because the Shannon-Nyquist criteria is not respected. Downsample() doesn't make sure the signal is filtered. Decimate does. So if need to choose, choose decimate ! +%sound(y_downsampled_ds,Fs_downsampled_ds) %Has distortion. This is because the Shannon-Nyquist criteria is not respected. Downsample() doesn't make sure the signal is filtered. Decimate does. So if need to choose, choose decimate ! end \ No newline at end of file