251 lines
8.9 KiB
Matlab
251 lines
8.9 KiB
Matlab
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 |