-% Cette routine permet de comparer deux méthodes de traitement
-% des sinus glissants (glissement logarithmique)
-% Le signal d'excitation est généré selon la méthode proposée par Farina
-% les deux méthodes de traitement sont :
-% - Fonction de Transfert
-% - Déconvolution par filtre inverse
-% Le système de test est un filtre passe-bas (type RC)
-% les paramètres à faire varier sont :
-% f1 et f2, les fréquences encadrant le signal généré
-% nwin le nombre de points sur lequel on fenêtre le signal (en début et en fin)
-% tmax : la durée du sinus glissant
-
-
-clear all
-close all
-%paramètres à faire varier
-f1=5; %fréquence de départ du sinus glissant (Hz)
-f2=8000; %fréquence d'arrivée du sinus glissant (Hz)
-nwin=400; % somme des nombres de point à fenêtrer en début et fin de signal
-tmax=0.5; %durée du sinus glissant (s)
-
-%Génération du chirp (Farina)
-fs=44100; %fréquence d'échantillonnage (Hz)
-dt=1/fs; %pas temporel
-t=0:dt:tmax; %vecteur temps
-om1=2*pi*f1; %pulsation basse
-om2=2*pi*f2; %pulsation haute
-K=tmax*om1/log(om2/om1);
-L=tmax/log(om2/om1);
-yc=sin(K*(exp(t./L)-1)); %Sinus glissant (Farina)
-
-%fenêtrage du chirp
-win=hanning(nwin);
-siz=length(yc);
-wint=ones(1,siz);
-wint(1:nwin/2)=win(1:nwin/2);
-wint(siz-nwin/2+1:siz)=win(nwin/2+1:nwin);
-y=yc.*wint; %Sinus glissant fenêtré
-% yinv=fliplr(y);
-% co=conv(y,yinv);
-% figure
-% plot(co)
-
-%test sur filtre passe-bas : contruction du signal filtré par méthode
-%récursive
-fc=400; %fréquence de coupuer du filtre
-RC=1/2/pi/fc;
-alf=dt/(dt+RC);
-yf(1)=y(1);
-for n=2:siz
- yf(n)=alf*y(n)+(1-alf)*yf(n-1); %Signal filtré
-end
-figure (1)
-plot(t,yf)
-xlabel('Time (s)')
-title('Input signal')
-grid on
-
-%Déconvolution par filtre inverse en temporel (Farina)
-yinv=fliplr(yc);
-fl=logspace(log10(f1),log10(f2),siz);
-yinv=yinv./fl;
-yinvf=yinv.*wint;
-IR=conv(yf,yinvf);
-figure(2)
-subplot(2,1,1)
-plot(IR)
-grid on
-xlabel('samples')
-title('Impulse response')
-ylabel('Deconvolution')
-sizi=length(IR);
-FR=fft(IR);
-FR=FR(1:sizi/2);
-fconv=0:fs/sizi:fs/2;
-fconv=fconv(1:sizi/2);
-LFR=20*log10(abs(FR));
-autoc=conv(yc,yinv);
-yeff=sqrt(sum(autoc.^2)*fs/2/(f2-f1));
-% yeff=sizi/2*(f2-f1)/fs
-% yeff=1
-sf=20*log10(yeff);
-LFR=LFR-sf;
-figure(3)
-subplot(2,1,1)
-% title('Frequency response (dB)')
-semilogx(fconv,LFR)
-title('Frequency response (dB)')
-axis([10 20000 -40 +5])
-grid on
-hold on
-FRT=abs(1./(1+j*2*pi*fconv*RC));
-LFRT=20*log10(FRT);
-semilogx(fconv,LFRT,'r')
-xlabel('Frequency (Hz)')
-legend('Processed','Theory',3)
-
-%Méthode de la fonction de transfert
-specy=fft(y);
-specyf=fft(yf);
-TF=specyf./specy;
-IR2=fftshift(ifft(TF));
-figure(2)
-subplot(2,1,2)
-plot(IR2)
-grid on
-xlabel('samples')
-ylabel('Inverse Transfer Function')
-TF2=TF(1:siz/2);
-f=0:fs/siz:fs/2;
-f=f(1:siz/2);
-LFR2=20*log10(abs(TF2));
-figure(3)
-subplot(2,1,2)
-semilogx(f,LFR2)
-axis([10 20000 -40 +5])
-grid on
-hold on
-FRT2=abs(1./(1+j*2*pi*f*RC));
-LFRT2=20*log10(FRT2);
-semilogx(f,LFRT2,'r')
-xlabel('Frequency (Hz)')
-legend('Processed','Theory',3)
+% Cette routine permet de comparer trois méthodes de traitement\r
+% des sinus glissants (glissement logarithmique)\r
+% Le signal d'excitation est généré selon la méthode proposée par Farina\r
+% les trois méthodes de traitement sont :\r
+% - Fonction de Transfert \r
+% - Déconvolution par filtre inverse\r
+% - Déconvolution par filtre inverse avec passage dans le domaine\r
+% fréquentiel\r
+% Le système de test est un filtre passe-bas (type RC)\r
+% les paramètres à faire varier sont :\r
+% f1 et f2, les fréquences encadrant le signal généré\r
+% nwin le nombre de points sur lequel on fenêtre le signal (en début et en fin)\r
+% tmax : la durée du sinus glissant\r
+\r
+clear all\r
+close all\r
+%input parameters\r
+f1=1; % start frequency (Hz)\r
+f2=10000; %end frequency(Hz)\r
+tmax=.5; % excitation signal duration (s)\r
+nwin=100; % windowing point number (half at the begining - half at the end)\r
+\r
+%chirp generation (Farina)\r
+fs=44100; % sampling frequency (Hz)\r
+dt=1/fs; % sampling period \r
+t=0:dt:tmax; % time vector\r
+om1=2*pi*f1; % begining pulsation\r
+om2=2*pi*f2; % ending pulsation \r
+K=tmax*om1/log(om2/om1);\r
+L=tmax/log(om2/om1);\r
+yc=sin(K*(exp(t./L)-1)); % excitation signal (logarithmic swept signal)\r
+\r
+%chirp windowing\r
+win=hanning(nwin); \r
+siz=length(yc);\r
+wint=ones(1,siz);\r
+wint(1:nwin/2)=win(1:nwin/2);\r
+wint(siz-nwin/2+1:siz)=win(nwin/2+1:nwin);\r
+y=yc.*wint; %windowed swept sine\r
+\r
+%Test on a low pass filter: processing of the filtered signal (iterative\r
+% method)\r
+fc=400; % cut-off frequency\r
+RC=1/2/pi/fc; \r
+alf=dt/(dt+RC);\r
+yf(1)=y(1);\r
+for n=2:siz\r
+ yf(n)=alf*y(n)+(1-alf)*yf(n-1); % filtered signal\r
+end\r
+\r
+%plot of the input and filtered signals\r
+figure (1)\r
+plot(t,y,'b'); \r
+xlabel('Time (s)');\r
+hold on \r
+plot(t,yf,'r')\r
+title('Input signals')\r
+grid on\r
+h = legend('Input signal','filtered signal',2);\r
+set(h,'location','SouthWest')\r
+\r
+% Impulse response: Inverse filter deconvolution (time domain)\r
+yinv=fliplr(y);\r
+fl=logspace(log10(f1),log10(f2),siz);\r
+yinv=yinv./fl; % inverse filter\r
+tic\r
+IR=conv(yf,yinv); % processed impulse response\r
+toc\r
+sizi=length(IR);\r
+is=fft(y).*conj(fft(yinv));\r
+norm=sqrt(sum(abs(is.^2))/siz);\r
+yeff=norm*sqrt(fs/2/(f2-f1));\r
+IR=IR*fs/yeff;\r
+t2=(-siz+1:1:siz-1).*dt;\r
+figure(2)\r
+subplot(3,1,1)\r
+plot(t2,IR,'b')\r
+grid on\r
+xlabel('samples')\r
+title('Impulse response')\r
+ylabel('Deconvolution')\r
+IRT=exp(-t2/RC)/RC;\r
+IRT(1:siz-1)=0;\r
+hold on \r
+plot(t2,IRT,'r')\r
+\r
+% frequency response : Inverse filter deconvolution (time domain)\r
+FR=fft(IR)/fs;\r
+FR=FR(1:sizi/2);\r
+fconv=0:fs/sizi:fs/2;\r
+fconv=fconv(1:sizi/2);\r
+LFR=20*log10(abs(FR));\r
+figure(3)\r
+subplot(3,1,1)\r
+semilogx(fconv,LFR)\r
+ylabel('Transfert function (decovonlution)')\r
+title('Frequency response (dB)')\r
+axis([10 20000 -40 +5])\r
+grid on\r
+hold on\r
+FRT=abs(1./(1+j*2*pi*fconv*RC));\r
+LFRT=20*log10(FRT);\r
+semilogx(fconv,LFRT,'r')\r
+xlabel('Frequency (Hz)')\r
+legend('Processed','Theory',3)\r
+\r
+%Transfer function method\r
+specy=fft(y);\r
+specyf=fft(yf);\r
+TF=specyf./specy;\r
+IR2=ifft(TF)*fs;\r
+figure(2)\r
+subplot(3,1,2)\r
+pad=zeros(1,siz-1);\r
+IR2=[pad IR2];\r
+plot(t2,IR2,'b')\r
+grid on\r
+hold on \r
+plot(t2,IRT,'r')\r
+xlabel('samples')\r
+ylabel('Inverse Transfer Function')\r
+TF2=TF(1:siz/2);\r
+f=0:fs/siz:fs/2;\r
+f=f(1:siz/2);\r
+LFR2=20*log10(abs(TF2));\r
+figure(3)\r
+subplot(3,1,2)\r
+semilogx(f,LFR2)\r
+axis([10 20000 -40 +5])\r
+grid on\r
+hold on\r
+FRT2=abs(1./(1+j*2*pi*f*RC));\r
+LFRT2=20*log10(FRT2);\r
+semilogx(f,LFRT2,'r')\r
+xlabel('Frequency (Hz)')\r
+ylabel('Transfer Function')\r
+legend('Processed','Theory',3)\r
+\r
+% Impulse response: Inverse filter deconvolution (freq domain)\r
+% yinv=fliplr(yc);\r
+% fl=logspace(log10(f1),log10(f2),siz);\r
+% yinv=yinv./fl; % inverse filter\r
+% yinvf=yinv.*wint; %windowed inverse filter\r
+% IR=conv(yf,yinvf); % processed impulse response\r
+tic \r
+iRF=fft([yf, zeros(1,siz-1)]).*fft([yinv, zeros(1,siz-1)])/yeff;\r
+IR3=real(ifft(iRF))*fs;\r
+toc\r
+% norm=sqrt(sum(abs(is.^2))/siz);\r
+% yeff=norm*sqrt(fs/2/(f2-f1));\r
+% scal=fs/yeff;\r
+t2=(-siz+1:1:siz-1).*dt;\r
+figure(2)\r
+subplot(3,1,3)\r
+plot(t2,IR3,'b')\r
+grid on\r
+xlabel('samples')\r
+title('Impulse response')\r
+ylabel('Deconvolution freq')\r
+IRT=exp(-t2/RC)/RC;\r
+IRT(1:siz-1)=0;\r
+hold on \r
+plot(t2,IRT,'r')\r
+\r
+% frequency response : Inverse filter deconvolution (freq domain)\r
+FR3=iRF(1:sizi/2);\r
+fconv=0:fs/sizi:fs/2;\r
+fconv=fconv(1:sizi/2);\r
+LFR3=20*log10(abs(FR3));\r
+figure(3)\r
+subplot(3,1,3)\r
+semilogx(fconv,LFR3)\r
+ylabel('Transfert function (deconv freq)')\r
+title('Frequency response (dB)')\r
+axis([10 20000 -40 +5])\r
+grid on\r
+hold on\r
+FRT=abs(1./(1+j*2*pi*fconv*RC));\r
+LFRT=20*log10(FRT);\r
+imin=round(f1/fs*sizi)+1;\r
+imax=round(f2/fs*sizi)+1;\r
+E3=sqrt(sum(abs(FRT(imin:imax)-FR3(imin:imax)).^2./(abs(FRT(imin:imax)).^2)))*100/(imax-imin+1);\r
+E3t = strcat('Erreur quadratique :',num2str(E3,2),' %');\r
+text(3000,-5,E3t);\r
+semilogx(fconv,LFRT,'r')\r
+xlabel('Frequency (Hz)')\r
+legend('Processed','Theory',3)\r
projects / cnaq.git / commitdiff