function speech_analysis()
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Task: To analyse the temporal and frequency plots of a signal.
%
% Inputs: -
%
% Outputs: -
%	
%
% Author: Charles Stelandre - charles.stelandre@ecam.fr
%
% Date: 09/04/2025
%
% Notes : This is the main file for the analysis of a signal. The main
% signal is modulator22.wav, which is present in the "sound" folder. Sound
% functions "sound()" are commented at then of the script.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all
close all
clc


% Construct the full file path
filename='modulator22.wav';
filepath = strcat('./sound/',filename);
% Read the audio file
[y, Fs] = audioread(filepath);
disp(['Successfully read the audio file: ', filename]);
disp(['Sampling frequency (Fs): ', num2str(Fs), ' Hz']);
disp(['Number of samples: ', num2str(length(y))]);

% Construct the output filename correctly
%[~, name, ~] = fileparts(filepath); % Get the filename without extension
%outputFilename = fullfile('.', ['processed_', name, '.wav']); % Create the new filename

% Write the audio to a new file with double the sampling rate
%audiowrite(outputFilename, y, Fs*2);
%disp(['Successfully wrote the processed audio to: ', outputFilename, ' with double the sampling rate.']);

disp('Playing the audio with double the sampling rate.');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Plot
t = (0:length(y)-1) / Fs; % Time in seconds

figure;
plot(t, y);
xlabel('Time (seconds)');
ylabel('Amplitude');
title(['Temporal Variation of ', filename]);
grid on;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Frequency Spectrum
%FFT

[yFFT, FFT_Time]=frequencySpectrum(y,Fs, 1);
disp("FFT duration :"); disp(FFT_Time);
%DFT

[yDFT, DFT_Time]=frequencySpectrum(y,Fs, 0);
disp("DFT duration :"); disp(DFT_Time);

%Modify the padding to make the change.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Spectrogram (Step Size = 5, Window Size = 50)
spectrogram(y, Fs, 5,50)
colorbar;
ylabel('Frequency (Hz)');
xlabel('Time (s)');
%% Spectrogram (Step Size = 30, Window Size = 50)
spectrogram(y, Fs, 30,50)
colorbar;
ylabel('Frequency (Hz)');
xlabel('Time (s)');
%% Spectrogram (Step Size = 5, Window Size = 5)
spectrogram(y, Fs, 5,5)
colorbar;
ylabel('Frequency (Hz)');
xlabel('Time (s)');
%% Spectrogram (Step Size = 30, Window Size = 5)
spectrogram(y, Fs, 1,50)
colorbar;
ylabel('Frequency (Hz)');
xlabel('Time (s)');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Going to Pratt, we see that :
%F0 : (100 + 130 +  100 + 120 + 100 + 90) / 6
%F1 : 578.3725189859462, 418.70239431349677, 552.8090680139439,
%308.88658136343446, 314.17710770594937, 363.8180262223959
%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() ---
downsample_factor_ds = round(Fs / desiredFreq); % Same result as = 6 for 4000Hz
%downsample_factor_ds=6;
y_downsampled_ds = downsample(y, downsample_factor_ds);
Fs_downsampled_ds = Fs / downsample_factor_ds;
disp(['--- Downsampling using downsample() ---']);
disp(['New sampling frequency (downsample): ', num2str(Fs_downsampled_ds), ' Hz']);
disp(['Number of samples (downsample): ', num2str(length(y_downsampled_ds))]);

% --- Downsampling using decimate() ---
downsample_factor_dec = round(Fs / desiredFreq); % Same result as = 6 for 4000Hz
%downsample_factor_dec=6;
y_decimated = decimate(y, downsample_factor_dec);
Fs_decimated = Fs / downsample_factor_dec;
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 ---
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);
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);
xlabel('Time (seconds)');
ylabel('Amplitude');
title(['Decimated Signal (decimate, Fs = ', num2str(Fs_decimated), ' Hz)']);
grid on;
%{
%% --- 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)
title('FFT of Decimated Signal (decimate)');
%}
%{
% --- Spectrograms of Downsampled Signals ---
figure;
subplot(2,1,1);
spectrogram(y_downsampled_ds, round(0.02*Fs_downsampled_ds), round(0.01*Fs_downsampled_ds), 512, Fs_downsampled_ds, 'yaxis');
title(['Spectrogram of Downsampled Signal (downsample, Fs = ', num2str(Fs_downsampled_ds), ' Hz)']);
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)']);
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
b_fir = fir1(order_fir, normalized_cutoff_fir, 'low'); % 'low' specifies a low-pass filter
a_fir = 1;
y_fir_filtered = filter(b_fir, a_fir, y); % Apply the FIR filter

figure;
freqz(b_fir,a_fir, 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');

% % FIR and 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
zeroPole(a_iir, b_iir,1);
zeroPole(a_fir, b_fir,1);

%% --- 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 !

end