// TP 4 /* PARTIE 1 : Implémentaion de la matrice de transformée de Fourier discrète */ // ------------ // Q1.(a) et Q1.(b) : voir le cours */ // ------------ /* Fonction permettant de débuter avec un indice =0 */ function j=I(i) j = i+1; endfunction // ------------ // Q2 Implémentation naîve de la TFD // ------------ // ------------ // Q2.(a) // function y=TFD1(N) omegaN = exp (%i*2*%pi/N); OMEGAN = zeros(N,N); for n=0:N-1 for k=0:N-1 OMEGAN(I(n),I(k)) =omegaN.^(n*k) end end y=conj(OMEGAN)/N; endfunction // ------------ // Q2.(b) //------------ /* Code non performant car deux boucles "for" imbriquées qui induisent une complexité algorithmique de l'ordre de N^2 */ //------------ /* Q3.(a) et Q3.(b) */ //------------ /*Méthode plus efficace pour implémenter la TFD à partir d'une matrice N*N obtenue par le produit matriciel d'un vecteur et de sa transposée. En effet, en définissant le vecteur w de taille N de k-ième composante (k=0,...,N-1) wk= sqrt(2pi/N) n exp(i pi/4) , le produit matriciel w*w^T fournit une matrice carrée N*N de terme général en (n,k) : (2pi)/N * nk exp(i pi /2) Remarquons que ((2pi)/N )* n*k *exp(i pi /2)= -i ((2pi)/N) * n* k et composons chacune des entrées de la matrice par la fonction exponentielle, la matrice résultante est alors Omega_n. */ //------------ /* Q3.(c) */ //------------ function y=TFD2(N) w= sqrt(2*%pi/N)*(0:N-1)*exp(%i*%pi/4).' z=w.' *w y=conj(exp(z))/N; endfunction //------------ /* Q4 */ //------------ N=100:100:400; tic() for i=1:length(N) TFD1(N(i)); end Res_TFD1=toc() tic(); for i=1:length(N) TFD2(N(i)); end Res_TFD2=toc() //------------ /* Q5 Il y a un bug ....cela reste trop lourd ... */ //------------ //------------------------------------------------ /* PARTIE 2 FFT et débruitage d''un signal */ //------------------------------------------------ //------------ /* Q6 */ //------------ T=1 function y= period1(x) f_0=50; f_1=120; y=sin(2*%pi*f_0.*x)+sin(2*%pi*f_1.*x); endfunction scf(1) N=1000; frequences=0:N;// abscisses des fréquences (N+1 fréquences)) x=frequences./N; // equi-espacés de 0.001 f=period1(x) // on observe le signal sur (N+1) points plot(x,f) title("fonction f périodique de période T=1") gcf().children.grid = color("grey70")*[1 1] // ajoute une grille //------------ /* Q7 TFD de la fonction "period1' et spectre d'amplitude */ //------------ L=1:floor(N/2);//on conserve seulement la moitie des frequences fhat=fft(f) // le fhat donne le vecteur Y*(N+1) (à un N+1 près) scf(2) plot(frequences(L),abs(fhat(L))/(N+1)) // plot(frequences,abs(fhat)/(N+1)) // sur toutes les fréquences title('spectre amplitude de fhat' ) gca().data_bounds(3:4) = [0 0.7]; // fixe le cadre des ordonnées gca().tight_limits = "on"; gcf().children.grid = color("grey70")*[1 1] // ajoute une grille //------------ /* Q8 Signal bruité */ //------------ funclean = f + 2.5*rand(x,"normal"); scf(3) subplot(2,1,1) plot(x,f) title('fonction f') gca().data_bounds(3:4) = [-10 5]; gca().tight_limits = "on"; gcf().children.grid = color("grey70")*[1 1] // ajoute une grille subplot(2,1,2) plot(x,funclean,'r') title('fonction funclean') gcf().children.grid = color("grey70")*[1 1] // ajoute une grille //------------ /* Q9 Spectre d'amplitude associé à f et à funclean */ //------------ frequences=(0:N); L=1:floor(N/2); fhat=fft(f) funclearhat=fft(funclean) scf(4) subplot(2,1,1) plot(frequences(L),abs(fhat(L))) title('spectre amplitude de fhat' ) gca().data_bounds(3:4) = [0 600]; gca().tight_limits = "on"; gcf().children.grid = color("grey70")*[1 1] // ajoute une grille subplot(2,1,2) plot(frequences(L),abs(funclearhat(L)),'r') title('spectre amplitude de funclearhat') gca().data_bounds(3:4) = [0 600]; gca().tight_limits = "on"; gcf().children.grid = color("grey70")*[1 1] // ajoute une grille //------------ /* Q10 Débruitage brut de funclean */ //------------ function y=debruitfunclearhat(r,ghat) N=length(ghat) L=1:floor(N/2) z=abs(ghat(L)) y=z.*(z>r) endfunction scf(5) subplot(3,1,1) plot(frequences(L),abs(fhat(L))) title('spectre amplitude de fhat' ) gca().data_bounds(3:4) = [0 600]; gca().tight_limits = "on"; gcf().children.grid = color("grey70")*[1 1] subplot(3,1,2) plot(frequences(L),abs(funclearhat(L)),'r') gca().data_bounds(3:4) = [0 600]; gca().tight_limits = "on"; title('spectre amplitude de funclearhat') gcf().children.grid = color("grey70")*[1 1] subplot(3,1,3) plot(frequences(L),debruitfunclearhat(250,funclearhat),'g') title('debruitage spectre amplitude de funclearhat') gca().data_bounds(3:4) = [0 600]; gca().tight_limits = "on"; gcf().children.grid = color("grey70")*[1 1] //------------ /* Q11 Filtre dans l'e domaine fréquentiel */ //------------ r=250 H =(abs(funclearhat)>r); fcleanhat=funclearhat.*H; //------------ /* Q12 Retour dans le domaine temporel */ //------------ fclean = ifft(fcleanhat); scf(6) subplot(3,1,1) plot(x,f) title('f' ) gca().data_bounds(3:4) = [-10 5]; gca().tight_limits = "on"; gcf().children.grid = color("grey70")*[1 1] subplot(3,1,2) plot(x, fclean,'r') title('fclean') gca().data_bounds(3:4) = [-10 5]; gca().tight_limits = "on"; gcf().children.grid = color("grey70")*[1 1] subplot(3,1,3) plot(x,funclean,'g') title('funclean') gca().data_bounds(3:4) = [-10 5]; gca().tight_limits = "on"; gcf().children.grid = color("grey70")*[1 1]