101 lines
2.4 KiB
Mathematica
101 lines
2.4 KiB
Mathematica
|
clear;
|
||
|
|
close all;
|
||
|
|
|
||
|
|
% returns sampling frequency in Hz and data
|
||
|
|
[y,Fs] = audioread('myson.wav');
|
||
|
|
|
||
|
|
% Fs = sampling frequency, 32000 for fluteircam.wav
|
||
|
|
lenW = 0.04*Fs; % window of lenW samples, i.e. 40 ms
|
||
|
|
|
||
|
|
df=0.9765625; %%% la dsp est calcul\'ee tous les df Hz
|
||
|
|
ff=-Fs/2:df:Fs/2; % length 32769 for fluteircam.wav
|
||
|
|
ffsub = ff(1:11:end); % compression x11, length 2979
|
||
|
|
|
||
|
|
tt = (0:length(y)-1)/Fs;
|
||
|
|
ttsub = (lenW/2:lenW:length(y))/Fs;
|
||
|
|
|
||
|
|
dsps = zeros(length(ffsub), floor((length(y)-lenW+1)/lenW));
|
||
|
|
|
||
|
|
% dsp fft
|
||
|
|
T = 1/Fs; % Sampling period
|
||
|
|
L = lenW; % Length of signal
|
||
|
|
t = (0:L-1)*T; % Time vector
|
||
|
|
fftsp = zeros(L/2+1, floor((length(y)-lenW+1)/lenW));
|
||
|
|
fftspCpx = zeros(L/2+1, floor((length(y)-lenW+1)/lenW));
|
||
|
|
f = Fs*(0:(L/2))/L;
|
||
|
|
|
||
|
|
for i = 0:floor((length(y)-lenW+1)/lenW)
|
||
|
|
|
||
|
|
% compute dsp AR
|
||
|
|
[~, ~, ~, ~, mydsp] = mylevinsondurbin(y(lenW*i+1:lenW*(i+1))',200,Fs);
|
||
|
|
dsps(:,i+1) = mydsp(1:11:end)'; % compression x11, length 2979
|
||
|
|
|
||
|
|
% compute dsp fft
|
||
|
|
myfft = fft(y(lenW*i+1:lenW*(i+1)));
|
||
|
|
P2 = abs(myfft/L);
|
||
|
|
P1 = P2(1:L/2+1);
|
||
|
|
P1(2:end-1) = 2*P1(2:end-1);
|
||
|
|
fftsp(:,i+1) = P1;
|
||
|
|
end
|
||
|
|
|
||
|
|
% take only positive frequencies for dsp
|
||
|
|
ffsubp = ffsub(1,(length(ffsub)-1)/2+1:end);
|
||
|
|
dspsp = dsps((length(dsps)-1)/2+1:end,:);
|
||
|
|
|
||
|
|
% plot
|
||
|
|
figure()
|
||
|
|
plot(tt,y)
|
||
|
|
xlabel('temps (s)')
|
||
|
|
ylabel('amplitude (u.a.)')
|
||
|
|
title('signal fluteircam')
|
||
|
|
|
||
|
|
% figure()
|
||
|
|
% surf(ttsub,ffsub,dsps,'EdgeColor','None');
|
||
|
|
% xlabel('temps (s)')
|
||
|
|
% ylabel('fréquences (Hz)')
|
||
|
|
% zlabel('amplitudes (u.a.)')
|
||
|
|
% title('Full DSP AR signal fluteircam')
|
||
|
|
|
||
|
|
figure()
|
||
|
|
surf(ttsub,ffsubp,dspsp,'EdgeColor','None');
|
||
|
|
xlabel('temps (s)')
|
||
|
|
ylabel('fréquences (Hz)')
|
||
|
|
zlabel('amplitudes (u.a.)')
|
||
|
|
title('DSP AR signal fluteircam')
|
||
|
|
|
||
|
|
figure()
|
||
|
|
imagesc(ttsub,ffsubp,dspsp)
|
||
|
|
xlabel('temps (s)')
|
||
|
|
ylabel('fréquences (Hz)')
|
||
|
|
title('Amplitude DSP AR signal fluteircam')
|
||
|
|
|
||
|
|
figure()
|
||
|
|
surf(ttsub,f,fftsp,'EdgeColor','None');
|
||
|
|
xlabel('temps (s)')
|
||
|
|
ylabel('fréquences (Hz)')
|
||
|
|
zlabel('amplitudes (u.a.)')
|
||
|
|
title('DSP FFT signal fluteircam')
|
||
|
|
|
||
|
|
figure()
|
||
|
|
imagesc(ttsub,f,fftsp)
|
||
|
|
xlabel('temps (s)')
|
||
|
|
ylabel('fréquences (Hz)')
|
||
|
|
title('Amplitude DSP FFT signal fluteircam')
|
||
|
|
|
||
|
|
|
||
|
|
% take max amplitude frequency
|
||
|
|
[maxDspsp, maxIndDspsp] = max(dspsp);
|
||
|
|
maxFfsubp = ffsubp(maxIndDspsp);
|
||
|
|
figure
|
||
|
|
plot(ttsub,maxFfsubp)
|
||
|
|
xlabel('temps (s)')
|
||
|
|
ylabel('fréquences (Hz)')
|
||
|
|
title('Frequency max DSP AR signal fluteircam')
|
||
|
|
|
||
|
|
[maxFftsp, maxIndFftsp] = max(fftsp);
|
||
|
|
maxF = f(maxIndFftsp);
|
||
|
|
figure
|
||
|
|
plot(ttsub,maxF)
|
||
|
|
xlabel('temps (s)')
|
||
|
|
ylabel('fréquences (Hz)')
|
||
|
|
title('Frequency max DSP FFT signal fluteircam')
|