This commit is contained in:
Gauthier BASSEREAU 2024-03-25 10:57:55 +00:00
parent bc5cc1a126
commit e9de737c0e
3 changed files with 131 additions and 10 deletions

View File

@ -1,4 +1,4 @@
function [power, duration] = frequencySpectrum(signal, fs, pad, plot)
function [power, duration] = frequencySpectrum(signal, fs, pad, verbose, rangefr)
%%%%%%%%%%%%%%%%%%
%function power = frequencySpectrum(signal, fs, pad)
%
@ -31,7 +31,10 @@ power = abs(y).^2/n; % power of the DFT
[val, ind] = max(power); % find the mx value of DFT and its index
if (plot)
rangemin = rangefr(1);
rangemax = rangefr(2);
if (verbose)
% plots
figure;
@ -53,9 +56,13 @@ if (plot)
plot(f,power, 'r');
xlabel('Frequency (Hz)')
ylabel('Power (a.u.)')
xlim([rangemin rangemax]);
subplot(1,3,3) % log frequency plot
plot(f,10*log10(power/power(ind)));
power_db = 10*log10(power/power(ind));
plot(f, power_db);
xlabel('Frequency (Hz)')
ylabel('Power (dB)')
xlim([rangemin rangemax]);
end

View File

@ -7,13 +7,18 @@ pkg load signal
#xlabel('Time(s)');
#ylabel('Signal Amplitude (normalized unit)');
#audiowrite("modifiedmodulator.wav",signal,Fs/2);
[signal, Fs] = audioread("modifiedmodulator.wav");
[signal, Fs] = audioread("Sound/modulator22.wav");
t=0:1/Fs:length(signal)/Fs - 1/Fs;
#figure; % Create a new figure
#plot(tt,signall);
#xlabel('Time(s)');
#ylabel('Signal Amplitude (normalized unit)');
#player=audioplayer(signal, Fs, 8)
#play(player);
##figure; % Create a new figure
##plot(t,signal);
##xlabel('Time(s)');
##ylabel('Signal Amplitude (normalized unit)');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Parameters for measurements
num_measurements = 100; % Number of measurements
@ -23,12 +28,12 @@ durations_fft = zeros(1, num_measurements);
for i = 1:num_measurements
% Measure time taken for DFT
tic;
[power_dft, duration_dft] = frequencySpectrum(signal, Fs, false, false);
[power_dft, duration_dft] = frequencySpectrum(signal, Fs, false, false, [100 2700]);
durations_dft(i) = duration_dft;
% Measure time taken for FFT with padding
tic;
[power_fft, duration_fft] = frequencySpectrum(signal, Fs, true, false);
[power_fft, duration_fft] = frequencySpectrum(signal, Fs, true, false, [100 2700]);
durations_fft(i) = duration_fft;
end
@ -45,6 +50,48 @@ fprintf('\n');
fprintf('Average duration for FFT (with padding): %f seconds\n', avg_duration_fft);
fprintf('Standard deviation for FFT (with padding): %f seconds\n', std_dev_fft);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Searching for formants without low pass filter
##start_time = 0.7; % start time in seconds
##end_time = 1.3; % end time in seconds
##
##start_index = round(start_time * Fs) + 1; % start index
##end_index = round(end_time * Fs) + 1; % end index
##
##cropped_signal = signal(start_index:end_index); % cropped signal
##
##[power_dft, duration_dft] = frequencySpectrum(cropped_signal, Fs, false, true, [0 160]);
#spectrogram(signal, Fs, 5, 30);
#spectrogram(signal, Fs, 5, 5);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%formants with low pass filter
N=8;
fc=2700;
[b,a]= butter(N,fc/(Fs/2));
freqz(b,a);
signal_filtered=filter(b,a,signal);
start_time = 0.75; % start time in seconds
end_time = 0.95; % end time in seconds
start_index = round(start_time * Fs) + 1; % start index
end_index = round(end_time * Fs) + 1; % end index
cropped_signal = signal_filtered(start_index:end_index); % cropped signal
[power_dft, duration_dft] = frequencySpectrum(cropped_signal, Fs, false, true, [0 2700]);

67
fir_iir.m Normal file
View File

@ -0,0 +1,67 @@
clc
clear all
close all
pkg load signal
samplingFreq = 500; % in Hz
N = samplingFreq; % number of samples for a 1s buffer of data
Nfir=7; % filter order
cutOffFrequency=10; % cut-off frequency
% sampling weighted window (hanning)
whann = window(@hanning,Nfir);
% Coefficients estimation using the inverse of the Fourier transform and truncation of the impulse response
n=0:1:(Nfir-1);
h=2*cutOffFrequency*sinc(2*(n-(Nfir-1)/2)*cutOffFrequency/samplingFreq)/samplingFreq; % Impulse response of the filter
hhann=(h.').*whann; % Weighting h by the Hanning window
hhann=hhann/sum(hhann); % Normalisation of the sum od the coeffs to be 1
figure
freqz(hhann,1,N,samplingFreq); % Frequency response of the impulse response weighted by the Hanning window
title('FIR Low-pass, order 6');
subplot(2,1,1)
ylim([-100 20])
[B, A] = butter(Nfir, cutOffFrequency/(samplingFreq/2));
figure;
freqz(B,A,N,samplingFreq)
title('IIR Low-pass, order 6');
subplot(2,1,1)
ylim([-100 20])
[Bc, Ac] = cheby1(Nfir, 10, cutOffFrequency/(samplingFreq/2));
figure;
freqz(Bc,Ac,N,samplingFreq)
title('IIR Low-pass, order 6');
subplot(2,1,1)
ylim([-100 20])
Z=roots(B) % Zeros of the transfer function
P=roots(A) % Poles of the transfer function
% Representation of zeros/poles in the complex plan
figure
zplane(Z,P) % Draw zeros and poles
title ('Zeros and poles of the transfer function of the IIR filter')
legend('zeros','poles')
grid on
% Representations of the impulse response of the IIR filter
figure
impz(B,A) % Draw the impulse response
title ('Impulse response of the IIR filter')
xlabel('Sample (n)')
grid on
% Representations of the impulse response of the FIR filter
figure
impz(hhann,1) % Draw the impulse response
title ('Impulse response of the FIR filter')
xlabel('Sample (n)')
grid on