回复 #15 yangzj 的帖子
不是呀,我的意思是我的图和songzy41的xk2a。jpg不一样,这个应该是表现了sinc的旁瓣了的,我的变换为db后不是这样的呢。我的变换坐标后是这样的
呵呵,我也在想这个问题.不知道songzy41老师是怎么做出来的.我做了下1024点采样,细化10倍的效果.
全局图
局部图
注意这里纵座标都没有取对数.
对数图
[ 本帖最后由 yangzj 于 2007-11-14 20:01 编辑 ]
回复 #18 yangzj 的帖子
不知你的细化是怎么作出来的,跟我的也不一样,呵呵,把程序传上来吧,谢谢啦 Nfft=Nmax ; fft2=fft(double(X11),Nfft); abs2=16*abs(fft2)/Nfft*2; %fft变换取点数 2的n次方把这一段改成
Nfft=Nmax;
n=0:Nfft-1;
k=0:0.1:Nfft/2; %频率间隔为0.1即细化十倍.
X11=double(X11);
fft2=(X11-mean(X11))*exp(-j*2*pi*n'*k/Nfft); %做DTFT
abs2=16*abs(fft2)/Nfft*2; %不知道为什么乘16
做图改成
f=(fs/Nfft)*k;
figure,plot(f,abs2);
[ 本帖最后由 yangzj 于 2007-11-14 21:27 编辑 ]
回复 #20 yangzj 的帖子
16其实我也不知道的,就是得出来的幅值跟原来的不一样,差了个16倍,这个在论坛里问过的,但没得到合理解释。 我对楼主在3楼的程序并没有作多大的修改,这里只给出求频谱那部分(前面部分都没有动)。因为有很大的直流分量,所以处理中先消除了直流分量。Nfft=1024; %fft变换取点数 2的n次方
subplot 211; plot(Xbar);
title('Xbar');
Xbmean=sum(Xbar)/Nfft;
Xbar=Xbar-Xbmean;
fft1=fft(,Nfft);
abs1=abs(fft1)/Nfft*2;
f=(fs/Nfft)*(0:(Nfft/2-1));
subplot 212;plot(f,abs1(1:512));
title('Spectrum');
ylabel('Intensity'); xlabel('Frequency');
figure
plot(f,20*log10(abs1(1:512))); grid;
title('Spectrum');
ylabel('Intensity(dB)'); xlabel('Frequency');
回复 #22 songzy41 的帖子
看到了去掉直流分量后的结果了,可不明白一个地方就是fft1=fft(,Nfft);,为什么要取xbar的91:1024个点然后采用补零的方式呢,我那个结果是在1:1024个点的条件下算出来的。 这语句:fft1=fft(,Nfft);
我并没有修改,在3楼的interferogramoerror.m程序中就是这样的,楼主提供的程序自已不清楚?
回复 #24 songzy41 的帖子
嘻嘻,看到了,我这个程序改了好几个版本,看错啦,谢谢啦。
页:
1
[2]