|
原帖由 wd963852 于 2009-1-10 10:45 发表
不太对 现在是频域的信号满足卷积的条件即:p_f*gr=target; 要求p_f ,我想用FFT的相关性质p_f=fft(ifft(target)/ifft(gr)),即 p_f=fft(x_t)=fft(z_t./y_t)求出P_f来 我前面描述的可能不太清楚 本来想画个框图的 ...
楼主在1楼的程序中,有几个错误,
1,gr和target从图上看出频率是从负值排到正值的,但在IFFT之前一定要把频率安排成从0~fs/2,然后再是负频率,即gr和target都要做fftshift后再进行IFFT变换。
2,gr在IFFT后y_t有些数值很小,有些数值还是由于计算误差造成的,它在分母上很容易引起更大的运算误差。
为此,我在程序中对z_t加了一个窗函数(如同FIR滤波器设计时窗函数法),这样就能尽量避免计算中引起更大的运算误差。程序如下:
fm=40;
fs=2*fm;
deltaf=1;
deltav=(-fs):deltaf:fs;
fs1=deltav;
gr=1./(1+((deltav-13.2)/3.3).^2);
L=length(fs1);
gr1=fftshift(gr);
fs2=0:deltaf:fs*2;
subplot(2,2,1);
plot(fs2,gr1,'r'); axis tight; grid;
title('gr1');
yy_t=ifft(gr1);
yytr=real(yy_t);
yyti=imag(yy_t);
t1=(0: L-1)/fs;
subplot(2,2,2);
plot(t1,yytr,'r',t1,yyti,'b'); axis tight; grid;
legend('real','imag');
title('yyt');
target=rectpuls(fs1,80);
target1=fftshift(target);
z_t=ifft(target1);
subplot(2,2,3);
plot(fs2,target1,'r'); axis tight; grid;
title('target1');
ztr=real(z_t);
zti=imag(z_t);
subplot(2,2,4)
plot(t1,ztr,'r',t1,zti,'b'); axis tight; grid;
legend('real','imag');
title('z_t');
figure(2);
win=hamming(41)';
zz_t=ztr;
yzt=zeros(1,L);
yzt(1:21)=zz_t(1:21).*win(21:41);
yzt(L-19:L)=zz_t(L-19:L).*win(1:20);
subplot(2,2,1);
plot(t1,yzt);axis tight; grid;
title('yzt');
x_t=zeros(1,L);
for k=1 : L
if yzt(k)~=0
x_t(k)=yzt(k)/yy_t(k);
end
end
xtr=real(x_t);
xti=imag(x_t);
subplot(2,2,2);
plot(t1,xtr,'r',t1,xti,'b');axis tight; grid;
legend('real','imag');
title('x_t');
p_f=fft(x_t);
p_f1=fftshift(p_f);
subplot(2,1,2); plot(fs1,abs(p_f1));
axis tight; grid;
title('p_f1'); |
|