diff --git a/rossignol/parametrique/fluteircam-dsp_AR.png b/rossignol/parametrique/fluteircam-dsp_AR.png new file mode 100644 index 0000000..03753bf Binary files /dev/null and b/rossignol/parametrique/fluteircam-dsp_AR.png differ diff --git a/rossignol/parametrique/fluteircam-dsp_AR_full.png b/rossignol/parametrique/fluteircam-dsp_AR_full.png new file mode 100644 index 0000000..999b9fd Binary files /dev/null and b/rossignol/parametrique/fluteircam-dsp_AR_full.png differ diff --git a/rossignol/parametrique/fluteircam-dsp_fft.png b/rossignol/parametrique/fluteircam-dsp_fft.png new file mode 100644 index 0000000..e55e269 Binary files /dev/null and b/rossignol/parametrique/fluteircam-dsp_fft.png differ diff --git a/rossignol/parametrique/fluteircam-dsp_max_AR.png b/rossignol/parametrique/fluteircam-dsp_max_AR.png new file mode 100644 index 0000000..ecdc3f6 Binary files /dev/null and b/rossignol/parametrique/fluteircam-dsp_max_AR.png differ diff --git a/rossignol/parametrique/fluteircam-dsp_max_fft.png b/rossignol/parametrique/fluteircam-dsp_max_fft.png new file mode 100644 index 0000000..65e7313 Binary files /dev/null and b/rossignol/parametrique/fluteircam-dsp_max_fft.png differ diff --git a/rossignol/parametrique/fluteircam-signal.png b/rossignol/parametrique/fluteircam-signal.png new file mode 100644 index 0000000..b67b84e Binary files /dev/null and b/rossignol/parametrique/fluteircam-signal.png differ diff --git a/rossignol/parametrique/mylevinsondurbin.m b/rossignol/parametrique/mylevinsondurbin.m new file mode 100644 index 0000000..1c4a670 --- /dev/null +++ b/rossignol/parametrique/mylevinsondurbin.m @@ -0,0 +1,57 @@ +%%% Algorithme de Levinson-Durbin pour la d\'etermination des param\`etres AR +%%% +%%% entr\'ees : +%%% - xx : signal +%%% - pp : ordre du mod\`ele AR (choisi de mani\`ere ind\'ependante) +%%% - fe : fr\'equence d'\'echantillonnage +%%% +%%% sorties : +%%% - aa : les param\`etres AR +%%% - sigma2 : variance du bruit +%%% - ref : coefficients de r\'eflexion +%%% - ff : fr\'equences auxquelles la dsp a \'et\'e calcul\'ee +%%% - mydsp : la dsp elle-m\^eme +%%% +%%% exemple : fe=32000;f0=440;xx=cos(2*pi*f0/fe*[1:1280]+2*pi*rand(1,1));mylevinsondurbin(xx,4,fe); +%%% +%%% S. Rossignol -- 2012 + +function [aa, sigma2, ref, ff, mydsp] = mylevinsondurbin (xx, pp, fe) + +acf = xcorr(xx, pp+1, 'biased'); %% autocorr\'elation +acf(1:pp+1) = []; %% on enl\`eve la partie n\'egative +acf(1) = real(acf(1)); %% Levinson-Durbin requiert c(1)==conj(c(1)) + +ref = zeros(pp,1); +gg = -acf(2)/acf(1); +aa = [ gg ]; +sigma2 = real( ( 1 - gg*conj(gg)) * acf(1) ); %% real : enl\`eve une \'eventuelle + %% partie imaginaire r\'esiduelle +ref(1) = gg; +for tt = 2 : pp + gg = -(acf(tt+1) + aa * acf(tt:-1:2)') / sigma2; + aa = [ aa + gg*conj(aa(tt-1:-1:1)), gg ]; + sigma2 = sigma2 * ( 1 - real(gg*conj(gg)) ); + ref(tt) = gg; +end; +aa = [1, aa]; + +%%% densit\'e spectrale de puissance +interm2=-j*2*pi/fe*[1:pp]; +df=0.9765625; %%% la dsp est calcul\'ee tous les df Hz +ff=-fe/2:df:fe/2; + +interm3=interm2'*ff; +interm=1.+aa(2:pp+1)*exp(interm3); +mydsp = sigma2./(interm.*conj(interm)); + +% figure(1); +% clf; +% grid on; +% hold on; +% plot(ff,mydsp,'linewidth',2); +% xlabel('frequency (in Hz)','fontsize',20); +% ylabel('magnitude','fontsize',20); +% hold off; +% drawnow; + diff --git a/rossignol/parametrique/myson-amp_dsp_AR.png b/rossignol/parametrique/myson-amp_dsp_AR.png new file mode 100644 index 0000000..7e1f8de Binary files /dev/null and b/rossignol/parametrique/myson-amp_dsp_AR.png differ diff --git a/rossignol/parametrique/myson-amp_dsp_AR_full.png b/rossignol/parametrique/myson-amp_dsp_AR_full.png new file mode 100644 index 0000000..5bc09c1 Binary files /dev/null and b/rossignol/parametrique/myson-amp_dsp_AR_full.png differ diff --git a/rossignol/parametrique/myson-amp_dsp_FFT.png b/rossignol/parametrique/myson-amp_dsp_FFT.png new file mode 100644 index 0000000..8f5d118 Binary files /dev/null and b/rossignol/parametrique/myson-amp_dsp_FFT.png differ diff --git a/rossignol/parametrique/myson-amp_dsp_FFT_full.png b/rossignol/parametrique/myson-amp_dsp_FFT_full.png new file mode 100644 index 0000000..5afdaf6 Binary files /dev/null and b/rossignol/parametrique/myson-amp_dsp_FFT_full.png differ diff --git a/rossignol/parametrique/myson-dsp_max_AR.png b/rossignol/parametrique/myson-dsp_max_AR.png new file mode 100644 index 0000000..29fb205 Binary files /dev/null and b/rossignol/parametrique/myson-dsp_max_AR.png differ diff --git a/rossignol/parametrique/myson-dsp_max_fft.png b/rossignol/parametrique/myson-dsp_max_fft.png new file mode 100644 index 0000000..f8322d6 Binary files /dev/null and b/rossignol/parametrique/myson-dsp_max_fft.png differ diff --git a/rossignol/parametrique/myson-signal.png b/rossignol/parametrique/myson-signal.png new file mode 100644 index 0000000..5a52c39 Binary files /dev/null and b/rossignol/parametrique/myson-signal.png differ diff --git a/rossignol/parametrique/parametrique.m b/rossignol/parametrique/parametrique.m new file mode 100644 index 0000000..7606651 --- /dev/null +++ b/rossignol/parametrique/parametrique.m @@ -0,0 +1,101 @@ +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') \ No newline at end of file