Merge branch 'son_param' into 'master'

Son param

See merge request hschindler/modelisationanalysespectrale!1

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;
+

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
+

rossignol/parametrique/parametrique_AR_fluteircam.m

@@ -0,0 +1,88 @@
+clear;
+close all;
+
+% returns sampling frequency in Hz and data
+[y,Fs] = audioread('fluteircam.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 = (1280/2:1280: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));
+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(blackman(lenW).*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()
+surf(ttsub,f,fftsp,'EdgeColor','None');
+xlabel('temps (s)')
+ylabel('fréquences (Hz)')
+zlabel('amplitudes (u.a.)')
+title('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')

rossignol/parametrique/parametrique_AR_myson.m

@@ -0,0 +1,71 @@
+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()
+imagesc(ttsub,ffsubp,dspsp)
+xlabel('temps (s)')
+ylabel('fréquences (Hz)')
+ylim([0 600])
+title('Amplitude DSP AR signal fluteircam')
+
+figure()
+imagesc(ttsub,f,fftsp)
+xlabel('temps (s)')
+ylabel('fréquences (Hz)')
+ylim([0 600])
+title('Amplitude DSP FFT signal fluteircam')

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');