求助:关于Matlab中Fourier变换问题.
大侠们好,我最近在使用Matlab编写反问题解法程序时遇到了麻烦,请大家不吝赐教!假设有这样的关系 F(w)=V(w)G(w).我这里用大写表示Fourier变换以后的函数。我已经知道了V的函数表达(假设为exp(-w)),g(t)的函数表达式也是知道的,这样就可以先对g(t)做Fourier变换然后再与V(w)做乘积即可。
但是Matlab是只能做离散Fourier变换的。而且我也只知道g(t)在的取值(t<0时可取0)。假设我在中取了100个点,做Fourier变换后的G(w)中w应该是0:2*pi:100*2*pi吧,然后我将V(w)也在同样的区间取值与G(w)做积。可是得到的却与我直接用f(t)做Fourier变换得出的结果不同(发现主要在前几个点上,第1个点相差太大,无法接受。。)。
不知道该如何解决此问题,请大侠们帮忙。 可能我的问题没有说清楚。我后面又在Matlab上演示发现有主要问题可能出在FFT算法上。
假设y=exp(-x^2);如果我先对y离散再进行离散Fourier变换(fft(nuy))和对y先Fourier变换(fourier(y))再离散结果相差太远,不知道是什么原因呢?
下面是我的程序:
clear
>> syms t y;
>> y=exp(-t^2);
>> tv=0:0.1:10;
>> nuy=subs(y,tv);
>> fft_nuy=fft(nuy);
>> fft_y=fourier(y);
>> xv=0:2*pi/100:2*pi;
>> plot(xv,abs(fft_nuy),'.')
>> hold on
>> nufft_y=subs(fft_y,xv);
>> plot(xv,abs(nufft_y),'r.-')
回复 2楼 的帖子
对fft我也有些不懂的地方,几个问题:
1. 我印象中,fft算法会周期延拓,所以不知道tv这样取可以否
2.fft算法后面一般要跟fftshift
3.fft后面好像有个啥系数,功率平方再除以数据长度之类的,我也没完全清楚
4.xv同tv要对应才有比较
把你程序稍做修改,并且把2条曲线归一化,发现是一模一样的。
clear
syms t y;
y=exp(-t^2);
tv=-10:0.1:10;
nuy=subs(y,tv);
fft_nuy=fftshift(fft(nuy));
fft_y=fourier(y);
N=length(tv);
xv=(-(N-1)/2:(N-1)/2)/N*1/0.1*2*pi;
% xv=0:2*pi/100:2*pi;
plot(xv,abs(fft_nuy)/max(abs(fft_nuy)),'*')
hold on
nufft_y=subs(fft_y,xv);
plot(xv,abs(nufft_y)/max(abs(nufft_y)),'r-')
回复 3楼 的帖子
非常感谢你的答复,依你的方法确实能得到两组相同的结果。但是我是想让tv的取值大于0时的结果。事实上我碰到是一个分段函数。我将t<0用0进行延拓。所以我对xv的取值很不确定(很抱歉我这方面的知识很欠缺,我从书上的理解是如果tv取值为0-1(fft似乎就是认为tv在其中取值),那么xv的取值是0:2*pi:(N-1)*2*pi)。
我把我实际要解决的问题的程序贴出来吧
%我要做的初步方案是检测g^=v1*f^,即g的Fourier变换后等于v1乘上f的Fourier变换
%这个等式已经被证明是正确的了。(当然上面的等式是针对连续F变换的)
syms xi t theta v1 v u g f
syms omega
x=0.30;
E=3;
tvalue=0.002:0.002:1;
xivalue=0:pi:500*2*pi;
%这个取值方法有待商榷,因为我在后面将tvalue从0-1延拓到了到了
%0-2所以我对xivalue的值是个什么样子不明白。。
delta=1.0e-3;
theta=sqrt(i*xi+0.25)-0.5;
u=(x+1)/t^1.5*exp(-(x+1-t)^2/(4*t));
g=2*sqrt(t)/t^2*exp(-(2-t)^2/(4*t)); %这个函数是我实际需要检测的函数
f=1/t^(1.5)*exp(-(1-t)^2/(4*t));
v=exp(-x*theta);
v1=exp(-theta);
nu_g=zeros(1,501);
nu_f=zeros(1,501);
nu_u=zeros(1,501);
nu_g(2:501)=subs(g,tvalue);
nu_f(2:501)=subs(f,tvalue);
nu_u(2:501)=subs(u,tvalue);
nu_v=subs(v,xivalue);
nu_v1=subs(v1,xivalue);
%-------------make periodization-------%
%下面的是做周期化处理也就是将0-1延拓到0-2使其
%具有周期性。这个方法我是参考国外一篇论文的方法。
spline_t_value=;
spline_g_value=;
spline_f_value=;
spline_u_value=;
nuspline_t=0:0.002:2;
pp=spline(spline_t_value,spline_g_value);
nuspline_g=;
pp=spline(spline_t_value,spline_f_value);
nuspline_f=;
pp=spline(spline_t_value,spline_u_value);
nuspline_u=;
normg=sqrt(intergral(,nu_g.^2));
%-------------make periodization-------%
fft_f=fft(nuspline_f);
fft_g=fft(nuspline_g);
cal_g=nu_v1.*fft_f;
%根据楼上的提示下面的也应该要修改
plot(xivalue,abs(fft_g),'r.-');
hold on
plot(xivalue,abs(cal_g),'-'); 不是很明白周期化延拓过程,这样做不是把原来的函数都变了么?
既然变了,那还有什么理由认为三者之间仍然满足原来的关系?
时域和频域的对应很简单的,
时域的采样点同频域的采样点个数相同,
时域两相邻采样点差的倒数为整个频域宽度,
同理,频域两相邻采样点差的倒数为整个时域宽度。
常用的代码如下:
ts=t(2)-t(1);
fs=1/ts;
f=(-N/2:N/2-1)/N*fs;
回复 5楼 的帖子
恩,应该是你所说的这样。我也很迷惑。但是文章上介绍的只是在t接近1时数据会受到影响。
页:
[1]