|
马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。
您需要 登录 才可以下载或查看,没有账号?我要加入
x
写了一个代码,自动生成信号
先采用时间域fft再进行空间域fft,其中时域fft采用了平均周期法,算出互谱,但出来的结果不太明白
在频域没有问题,
问题可能是,1.在第一次fft时是否需要取abs然后再进行第二次fft?
2.空间上的间距与空间频率的关系不太清楚?
%auto generate signal
sf=10000;
t=0:1/sf:1-1/sf;
dx=0.025;
x=0:dx:0.5;
f=300;
omega=2*pi*f;
c=340;
for j=1:length(x)
A(:,j)=10*sin(omega.*(t-x(j)/c));
end
Pressuref=[];
for i=2:21
%loading pressure
micone=A(:,1);
micother=A(:,i);
%the first fft for time domain
Ns= 1024;
NNs=Ns/2;
n=[0: Ns/10-1];
f = sf*n/Ns;
Nn=ceil((length(micone)-Ns)/NNs);
pxx2=0;
for in=0:Nn;
pxx2=pxx2+abs(fft(micone(in*NNs+1:min(Ns+NNs*in,length(micone))), Ns).*conj(fft(micother(in*NNs+1:min(Ns+NNs*in,length(micone))),Ns)))/Ns;
end
Pxx= 10 * log10((pxx2)/(in+1));
Pressuref=[Pressuref Pxx]; %the frequency matrix after fft
plot(f,Pxx(n+1)); grid;%plot the power spectrum for validate
end
%the second fft for spatial domain
PressureFK=[];
for jf=1:length(f)
PressureF=Pressuref(jf,:);
nfft=512;
PressureK=fft(PressureF,nfft);
mag=abs(PressureK)*2/nfft;
wavenumber=(0:(length(PressureK)/2-1))*dx/length(PressureK);
PressureFK=[PressureFK; mag(1:nfft/2)];
end
[X,Y]=meshgrid(f,wavenumber);
plot3(X,Y,PressureFK');grid ;
高人指点下吧
[ 本帖最后由 eight 于 2008-5-6 18:15 编辑 ] |
-
|