Fixed filters. Created plots for filtered compared to original outputs.

iirFilter.m

@@ -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;
-

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