diff --git a/rossignol/parametrique/fluteircam-dsp_max_MUSIC_p4M210.svg b/rossignol/parametrique/fluteircam-dsp_max_MUSIC_p4M210.svg new file mode 100644 index 0000000..e70183b --- /dev/null +++ b/rossignol/parametrique/fluteircam-dsp_max_MUSIC_p4M210.svg @@ -0,0 +1,216 @@ + + +02468101214temps (s)020040060080010001200140016001800fréquences (Hz)Frequency max DSP MUSIC signal flureircam (p = 4, M = 210) diff --git a/rossignol/parametrique/fluteircam-spectre_max_fft.png b/rossignol/parametrique/fluteircam-spectre_max_fft.png new file mode 100644 index 0000000..f9287f7 Binary files /dev/null and b/rossignol/parametrique/fluteircam-spectre_max_fft.png differ diff --git a/rossignol/parametrique/mymusic_matlab.m b/rossignol/parametrique/mymusic_matlab.m new file mode 100644 index 0000000..40afa07 --- /dev/null +++ b/rossignol/parametrique/mymusic_matlab.m @@ -0,0 +1,75 @@ +%%% Algorithme de Music pour la d\'etermination des param\`etres Music +%%% +%%% entr\'ees : +%%% - xx : signal +%%% - pp : ordre du mod\`ele (choisi de mani\`ere ind\'ependante) +%%% - MM : nombre de coefficients de corr\'elation pris en compte (MM>=pp) +%%% - fe : fr\'equence d'\'echantillonnage +%%% +%%% sorties : +%%% - ff : fr\'equences auxquelles la dsp a \'et\'e calcul\'ee +%%% - mydsp : la dsp elle-m\^eme +%%% +%%% exemples : clear;rand('seed',100*sum(clock));fe=32000;f0=440;tsig=1280;xx=cos(2*pi*f0/fe*[1:tsig]+2*pi*rand(1,1));mymusic(xx,2,10,fe); +%%% clear;rand('seed',100*sum(clock));fe=32000;f0=440;tsig=1280;xx=cos(2*pi*f0/fe*[1:tsig]+2*pi*rand(1,1))+cos(2*pi*(f0+26)/fe*[1:tsig]+2*pi*rand(1,1));mymusic(xx,4,300,fe); +%%% +%%% S. Rossignol -- 2012 + + +%%% utilisation en script : +%clear;rand('seed',100*sum(clock));fe=32000;f0=440;tsig=1280;xx=cos(2*pi*f0/fe*[1:tsig]+2*pi*rand(1,1))+cos(2*pi*(f0+26)/fe*[1:tsig]+2*pi*rand(1,1));pp=4;MM=400; +%clear;rand('seed',100*sum(clock));fe=32000;f0=440;tsig=1280;xx=cos(2*pi*f0/fe*[1:tsig]+2*pi*rand(1,1))+1e-2*randn(1,tsig);pp=2;MM=10; + + + +function [ff, mydsp] = mymusic_matlab(xx, pp, MM, fe) + + +if (MM<=pp) + fprintf(1, 'Il faut absolument MM>pp !!!\n'); + return; +end; + +MM1=MM-1; +res=xx; +xx = xx-mean(xx); + + +%%% corr\'elations +acf = xcorr(xx, MM1, 'biased'); +lMM=length(acf); +rrr1 = acf(MM1+1:lMM)'; +for ii=1:MM1 + rrr1 = [rrr1 acf(MM1+1-ii:lMM-ii)']; +end; +rrr1 = rrr1'; + + +%%% m\'ethode directe pour trouver toutes les valeurs propres +[v, lambda] = eig(rrr1); +lamb = diag(lambda); +[vl,pl] = sort(abs(lamb),'descend'); + + +%%% densit\'e spectrale de puissance +df=0.9765625; %%% la dsp est calcul\'ee tous les df Hz +ff=-fe/2:df:fe/2; + +mydsp=zeros(length(ff),1); +deni=zeros(MM,MM); +for ii=pp+1:MM + deni = deni + v(:,pl(ii))*conj(v(:,pl(ii)))'; +end; +for ii=1:length(ff) + ee = cos(2*pi*ff(ii)*[0:MM1]/fe); + den = conj(ee)*deni*ee'; + mydsp(ii) = abs(1/den); +end; + + +%%% on enl\`eve \'eventuellement une composante non nulle +mydsp = mydsp-min(mydsp); +%mydsp = mydsp/max(mydsp); %%% si on fait \c{c}a, c'est norme 1 + +end + diff --git a/rossignol/parametrique/myson-fft.png b/rossignol/parametrique/myson-fft.png new file mode 100644 index 0000000..bd9e54d Binary files /dev/null and b/rossignol/parametrique/myson-fft.png differ diff --git a/rossignol/parametrique/myson-fft_zoom1.png b/rossignol/parametrique/myson-fft_zoom1.png new file mode 100644 index 0000000..70e2541 Binary files /dev/null and b/rossignol/parametrique/myson-fft_zoom1.png differ diff --git a/rossignol/parametrique/myson-fft_zoom2.png b/rossignol/parametrique/myson-fft_zoom2.png new file mode 100644 index 0000000..680839c Binary files /dev/null and b/rossignol/parametrique/myson-fft_zoom2.png differ diff --git a/rossignol/parametrique/parametrique_fluteircam.m b/rossignol/parametrique/parametrique_AR_fluteircam.m similarity index 100% rename from rossignol/parametrique/parametrique_fluteircam.m rename to rossignol/parametrique/parametrique_AR_fluteircam.m diff --git a/rossignol/parametrique/parametrique.m b/rossignol/parametrique/parametrique_AR_myson.m similarity index 100% rename from rossignol/parametrique/parametrique.m rename to rossignol/parametrique/parametrique_AR_myson.m diff --git a/rossignol/parametrique/parametrique_MUSIC_fluteircam.m b/rossignol/parametrique/parametrique_MUSIC_fluteircam.m new file mode 100644 index 0000000..8041e50 --- /dev/null +++ b/rossignol/parametrique/parametrique_MUSIC_fluteircam.m @@ -0,0 +1,82 @@ + + +%rand('seed',100*sum(clock)); +%fe=32000; +%f0=440; +%tsig=1280; +%xx=cos(2*pi*f0/fe*[1:tsig]+2*pi*rand(1,1))+cos(2*pi*(f0+40)/fe*[1:tsig]+2*pi*rand(1,1))+cos(2*pi*2*f0/fe*[1:tsig]+2*pi*rand(1,1))+cos(2*pi*(2*f0+40)/fe*[1:tsig]+2*pi*rand(1,1)); +%xx=cos(2*pi*f0/fe*[1:tsig]+2*pi*rand(1,1))+1e-2*randn(1,tsig); +%pp=2; +%MM=10; + +% returns sampling frequency in Hz and data +[xx,fe] = audioread('fluteircam.wav'); +xx = xx'; + +pp=4; +MM=210; + +lenW = 0.04*fe; + +maxFreqs = []; + +for i = 0:floor((length(xx)-lenW+1)/lenW) + xxsub = xx(1,lenW*i+1:lenW*(i+1)); + [ff, mydsp] = mymusic_matlab(xxsub,pp,MM,fe); + + mydspP = mydsp(floor(length(mydsp)/2):end,:); + ffP = ff(:,floor(length(mydsp)/2):end); + [maxi,ind] = max(mydspP); + maxFreq = ffP(ind); + maxFreqs = [maxFreqs maxFreq]; +end + +figure() +tt = (lenW/2:lenW:length(xx))/fe; +plot(tt,maxFreqs) +xlabel('temps (s)') +ylabel('fréquences (Hz)') +title('Frequency max DSP MUSIC signal flureircam (p = 2, M = 200)') + +%%% figures +% figure(); +% clf; +% grid on; +% hold on; +% plot(ff,mydsp,'linewidth',2); +% xlabel('frequency (in Hz)','fontsize',20); +% ylabel('magnitude','fontsize',20); +% title('zoom MUSIC'); +% xlim([400 506]); +% hold off; + +% figure(); +% clf; +% grid on; +% hold on; +% plot(ff,mydsp,'linewidth',2); +% xlabel('frequency (in Hz)','fontsize',20); +% ylabel('magnitude','fontsize',20); +% title('MUSIC'); +% hold off; +% drawnow; + +% % montrer intéret paramétrique par rapport à FFT +% figure(); +% fftxx = abs(fftshift(fft(xx))); % fftshift permet de passer de la fft entre 0 et fe à la fft entre -fe/2 et fe/2 +% freq = linspace(-fe/2,fe/2,length(fftxx)); +% plot(freq,fftxx); +% grid on; +% xlabel('fréqunce (Hz)'); +% ylabel('amplitude (u.a.)'); +% title('fft'); +% +% figure(); +% xx = [xx zeros(1,32768-length(xx))]; % 2^15 = 32768 +% fftxx = abs(fftshift(fft(xx))); +% freq = linspace(-fe/2,fe/2,length(fftxx)); +% plot(freq,fftxx); +% grid on; +% xlabel('fréqunce (Hz)'); +% ylabel('amplitude (u.a.)'); +% title('fft zero padding'); \ No newline at end of file