Comment filtrer un cepstre

[Le_ramo]

Bonjour ,
j'ai un signal additionné à une version décalée de lui-même ,et je vais enlever cet écho pour obtenir le signal initial.
En effet j'ai utilisé l'analyse spectrale dans ce but.
j'ai un spectre qu'il faut filtrer pour éliminer les impulsions de celui-ci pour pouvoir reconstruire le signal original après.

Mon code est celui-ci :

Code :

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
% Test de l'analyse cepstrale sur les signaux déterministes
% Créé le 12/26/2011
% -------------------------------------------------------------------------

clear all;
clear data;
close all;
clc;

% 1- Définition des paramètres et du signal à étudier
%    ------------------------------
f1=50;
f2=150; 
N=1024; % Nbre de points du signal
tmin=-1;
tmax=1;
t=linspace(tmin,tmax,N); % Définition du l'axe de temps
fe=N/2;
fedemi=fe/2;
f=linspace(-fedemi,fedemi,N); % Définition du l'axe de fréquence
x1=cos(2*pi*f1*t);
x2=sin(2*pi*f2*t);
x_non_decale=x1.*x2;

%1.1- Signal normal
%     -----------------
figure;plot(t,x_non_decale);  % figure 1
xlabel('Temps t');
ylabel('x(t)');
legend('Signal original: x(t)=cos(2*pi*f1*t)*sin(2*pi*f2*t)');


%1.2- Signal avec écho 
%     -----------------
a0=0.5;
t0=0.2;
x_decale=cos(2*pi*f1*(t-t0)).*sin(2*pi*f2*(t-t0));

x_total=x_non_decale+a0*x_decale;
figure;plot(t,x_total);  % figure 2
xlabel('Temps t');
ylabel('x(t)+a0*x(t-t0)');
legend('Version décalée du signal x(t)');

% -------------------------------------------------------------------------

% 2- Application de l'analyse cepstrale sur le signal

%xf_non_decale=fftshift(fft(x_non_decale));
%xf_total2=xf_non_decale.*(1+a0*exp(-2*pi*j*f*t0));
xf_total=fftshift(fft(x_total));
figure;plot(f,xf_total);  %figure 3

module_xf_au_carre=abs(xf_total).^2;
%module_xf_au_carre2=(abs(xf_non_decale).^2).*(1+a0^2+2*a0*cos(2*pi*f*t0));

figure;plot(f,module_xf_au_carre);  %figure 4
%figure;plot(f,module_xf_au_carre2);  

log_module_xf_au_carre=log(module_xf_au_carre);
figure;plot(f,log_module_xf_au_carre);  %figure 5

cepstre_x=real(ifftshift(ifft(log(xf_total))));
figure;plot(t,cepstre_x,'b',t,real(ifftshift(ifft(log(1+a0^2+2*a0*cos(2*pi*f*t0))))),'r');   %figure 6

En effet, je veux enlever les deux impulsions rouges qui apparaissent dans la figure 6 pour pouvoire reconstruire le signal original à partir du cepstre liftré.

Comment choisir ma fenêtre de liftrage ??

[Le_ramo]

Nouveau code :

Code :

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
% Test de l'analyse cepstrale sur les signaux aléatoires
% Créé le 12/26/2011
% -------------------------------------------------------------------------

clear all;
clear data;
close all;
clc;

% 1- Définition des paramètres et du signal à étudier
%    ------------------------------
f1=50;
f2=150; 
N=1024;
tmin=0;
tmax=1;
t=linspace(tmin,tmax,N);
fe=N/2;
fedemi=fe/2;
f=linspace(0,fedemi,N);
x1=cos(2*pi*f1*t);
x2=sin(2*pi*f2*t);
x_non_decale=x1.*x2 +0.1*randn(size(t));

%1.1- Signal normal avec bruit
%     ------------------------
figure;plot(t,x_non_decale);  % figure 1
xlabel('Temps t');
ylabel('x(t)');
legend('Signal original: x(t)=cos(2*pi*f1*t)*sin(2*pi*f2*t)+bruit');


%1.2- Signal avec écho 
%     -----------------
a0=0.5;
t0=0.2;
x_decale=cos(2*pi*f1*(t-t0)).*sin(2*pi*f2*(t-t0));

x_total=x_non_decale+a0*x_decale;
figure;plot(t,x_total);  % figure 2
xlabel('Temps t');
ylabel('x(t)+a0*x(t-t0)');
legend('Version décalée du signal x(t)');

% -------------------------------------------------------------------------


% 2- Application de l'analyse cepstrale sur le signal
%    ------------------------------------------------


% 2.1- Application de l'algorithme  par la 1ère méthode (cepstre du signal)
%      --------------------------------------------------------------------

%xf_total=real(fftshift(fft(x_total)));
xf_total=real(fftshift(fft(x_total)));
figure;plot(f,xf_total);  %figure 3
xlabel('Fréquence f');
ylabel('xtotal(f)');
legend('Spectre du signal total');


module_xf_au_carre=abs(xf_total).^2;
%module_xf_au_carre2=(abs(xf_non_decale).^2).*(1+a0^2+2*a0*cos(2*pi*f*t0));

figure;plot(f,module_xf_au_carre);  %figure 4
xlabel('Fréquence f');
ylabel('xtotal(f)^2');
legend('Spectre du module au carré du signal total');

 

log_module_xf_au_carre=log(module_xf_au_carre);
figure;plot(f,log_module_xf_au_carre);  %figure 5
xlabel('Fréquence f');
ylabel('log_module_xf_au_carre');
legend('Cepstre du signal total');


cepstre_x=real(ifftshift(ifft(log(xf_total))));
figure;plot(t,cepstre_x,'b',t,real(ifftshift(ifft(log(1+a0^2+2*a0*cos(2*pi*f*t0))))),'r');   %figure 6
xlabel('Temps t');
ylabel('cepstre_x(t) et Cepstre_écho(t)');
legend('Cepstre du signal total','Cepstre des impulsions du écho');



bruit=real(ifftshift(ifft(log(1+a0^2+2*a0*cos(2*pi*f*t0)))));
cepstre_x_total=cepstre_x + bruit;

% 2.1.1- Choix du liftre 
%      ---------------
filtre=ones(1,1024);

filtre(1,200)=0;
filtre(1,306:308)=0;
filtre(1,408:412)=0;
filtre(1,612:618)=0;
filtre(1,716:719)=0;
filtre(1,816)=0;

% 2.1.2- liftrage du signal 
%      ------------------
cepstre_filtre=cepstre_x_total.*filtre;

figure;plot(t,cepstre_x);   %figure 7
xlabel('Temps t');
ylabel('cepstre_x(t)');
legend('Cepstre du signal sans écho');


figure;plot(t,cepstre_filtre);  %figure 8
xlabel('Temps t');
ylabel('cepstre_x(t)');
legend('Cepstre du signal filtré');


signal_filtre_f=real(fftshift(fft(exp(cepstre_filtre))));
signal_filtre=real(ifftshift(ifft(signal_filtre_f)));
figure;plot(t,x_non_decale,'b');  %figure 9
xlabel('Temps t');
ylabel('x(t)');
legend(' signal x(t) initial sans écho');

figure;plot(t,signal_filtre,'r');  %figure 10
xlabel('Temps t');
ylabel('xfiltré(t)');
legend(' signal x(t) après application du filtre sans écho');

axis([0 1 -1.2 1.2]);

Mais je n'arrive pas à retrouver le signal initial exact
Ceci à cause des paramètres du filtre que je n'arrive pas à préciser exactement !
Pas de solutions chers amis ?

[Winjerome]

Bonjour,

• Attention d'utiliser la fonction IFFTSHIFT avant IFFT dans les lignes (j'en oublie peut-être)
  
  Code :

  1
  2
  3
  cepstre_x=real(ifftshift(ifft(log(xf_total))));
  bruit=real(ifftshift(ifft(log(1+a0^2+2*a0*cos(2*pi*f*t0)))));
  signal_filtre=real(ifftshift(ifft(signal_filtre_f)));

• Pour la visualisation de la FFT du signal, on utilise plutôt la valeur absolue:
  
  Code :

  1
  2
  xf_total=real(fftshift(fft(x_total)));
  figure;plot(f,xf_total);  %figure 3

  devient
  
  Code :

  1
  2
  xf_total=fftshift(fft(x_total));
  figure;plot(f,abs(xf_total));  %figure 3

• Les fréquences à la sortie de la FFT sont situées dans [-fe/2 fe/2] donc
  
  Code :

  1
  2
  f=linspace(0,fedemi,N);  % dans [0 fe/2]
  bruit=real(ifftshift(ifft(log(1+a0^2+2*a0*cos(2*pi*f*t0)))));

  est incorrect.
• Code :

  signal_filtre_f=real(fftshift(fft(exp(cepstre_filtre))));

  La fonction EXP est à appliquer après la FFT.
• Tu prends un peu trop souvent la partie réelle, ce qui supprime une partie de l'information.

[Le_ramo]

Merci pour tes conseils.

Comment je peux filter mon cepstre des impulsions rouges qu'on voit sur la figure 6?
J'ai choisis un filtre =ones(1,1024) qu'il faut forcer à zéros ses éléments qui correspondent aux positions des impulsions dans la figure 6.
J'ai fais cela par main, c'est à dire j'ai utilisé la souris et le bouton "data cursor" pour capter les positions de ces impulsions, mais inutile, toujours pas de réponse exacte !

Prière de faire copier/coller du code et l'éxécuter dans Matlab pour voir la figure 6 et les impulsions...

Cordialement

[Winjerome]

Envoyé par Le_ramo

Prière de faire copier/coller du code et l'éxécuter dans Matlab pour voir la figure 6 et les impulsions...

Prière de fournir ton nouveau code avec les erreurs citées corrigées