This commit is contained in:
Gauthier BASSEREAU 2024-03-26 10:16:15 +00:00
parent c32df544bb
commit b8e8de1a9a
2 changed files with 85 additions and 33 deletions

View File

@ -58,6 +58,7 @@ if (verbose)
ylabel('Power (a.u.)')
xlim([rangemin rangemax]);
subplot(1,3,3) % log frequency plot
power_db = 10*log10(power/power(ind));
plot(f, power_db);
@ -65,4 +66,51 @@ if (verbose)
ylabel('Power (dB)')
xlim([rangemin rangemax]);
subplot(1,3,3) % log frequency plot
power_db = 10*log10(power/power(ind));
plot(f, power_db);
xlabel('Frequency (Hz)')
ylabel('Power (dB)')
xlim([rangemin rangemax]);
% Shift power_db to all positive values
shift_amount = min(power_db) * -1;
power_db_shifted = power_db + shift_amount;
% Find peaks and their indices
[pks, locs] = findpeaks(power_db_shifted);
% Interpolate to get a smooth line
xi = linspace(min(f(locs)), max(f(locs)), 1000);
yi = interp1(f(locs), pks, xi, 'pchip');
% Shift the interpolated values back down and up by 10
yi = yi - shift_amount + 10;
% Smooth the line using a moving average filter
windowSize = 15; % Adjust this value to change the amount of smoothing
b = (1/windowSize)*ones(1,windowSize);
yi = filter(b, 1, yi);
yi = filter(b, 1, flip(yi)); % Apply the filter in reverse direction
yi = flip(yi); % Flip the data back to original direction
% Shift yi to all positive values
shift_amount_yi = min(yi) * -1;
yi_shifted = yi + shift_amount_yi;
% Find peaks of the smoothed line
[pks_smooth, locs_smooth] = findpeaks(yi_shifted);
% Plot the envelope
hold on;
plot(xi, yi, 'r-'); % plot the envelope
% Plot points and labels at the peaks
for i = 1:length(pks_smooth)
plot(xi(locs_smooth(i)), pks_smooth(i) - shift_amount_yi, 'ko');
text(xi(locs_smooth(i)), pks_smooth(i) - shift_amount_yi, ['F' num2str(i)], 'VerticalAlignment', 'bottom');
end
hold off;
end

View File

@ -20,35 +20,35 @@ t=0:1/Fs:length(signal)/Fs - 1/Fs;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Parameters for measurements
num_measurements = 100; % Number of measurements
durations_dft = zeros(1, num_measurements);
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, [100 2700]);
durations_dft(i) = duration_dft;
% Measure time taken for FFT with padding
tic;
[power_fft, duration_fft] = frequencySpectrum(signal, Fs, true, false, [100 2700]);
durations_fft(i) = duration_fft;
end
% Calculate average and standard deviation
avg_duration_dft = mean(durations_dft);
std_dev_dft = std(durations_dft);
avg_duration_fft = mean(durations_fft);
std_dev_fft = std(durations_fft);
fprintf('Average duration for DFT: %f seconds\n', avg_duration_dft);
fprintf('Standard deviation for DFT: %f seconds\n', std_dev_dft);
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);
##% Parameters for measurements
##num_measurements = 100; % Number of measurements
##durations_dft = zeros(1, num_measurements);
##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, [100 2700]);
## durations_dft(i) = duration_dft;
##
## % Measure time taken for FFT with padding
## tic;
## [power_fft, duration_fft] = frequencySpectrum(signal, Fs, true, false, [100 2700]);
## durations_fft(i) = duration_fft;
##end
##
##% Calculate average and standard deviation
##avg_duration_dft = mean(durations_dft);
##std_dev_dft = std(durations_dft);
##
##avg_duration_fft = mean(durations_fft);
##std_dev_fft = std(durations_fft);
##
##fprintf('Average duration for DFT: %f seconds\n', avg_duration_dft);
##fprintf('Standard deviation for DFT: %f seconds\n', std_dev_dft);
##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
@ -72,19 +72,23 @@ fprintf('Standard deviation for FFT (with padding): %f seconds\n', std_dev_fft);
N=8;
fc=2700;
[b,a]= butter(N,fc/(Fs/2));
freqz(b,a);
%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
%One : 0.75-0.95 / 0.7-1.2
%Two : 1.3-1.7
%Three : 1.85-2.2
start_time = 2.0; % start time in seconds
end_time = 2.2; % 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]);
[power_dft, duration_dft] = frequencySpectrum(cropped_signal, Fs, false, true, [100 3000]);