声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 2791|回复: 5

[图像处理] 求助:关于Matlab中Fourier变换问题.

[复制链接]
发表于 2008-6-27 15:52 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有账号?我要加入

x
大侠们好,我最近在使用Matlab编写反问题解法程序时遇到了麻烦,请大家不吝赐教!
  假设有这样的关系 F(w)=V(w)G(w).我这里用大写表示Fourier变换以后的函数。我已经知道了V的函数表达(假设为exp(-w)),g(t)的函数表达式也是知道的,这样就可以先对g(t)做Fourier变换然后再与V(w)做乘积即可。
  但是Matlab是只能做离散Fourier变换的。而且我也只知道g(t)在[0,1]的取值(t<0时可取0)。假设我在[0,1]中取了100个点,做Fourier变换后的G(w)中w应该是0:2*pi:100*2*pi吧,然后我将V(w)也在同样的区间取值与G(w)做积。可是得到的却与我直接用f(t)做Fourier变换得出的结果不同(发现主要在前几个点上,第1个点相差太大,无法接受。。)。
  不知道该如何解决此问题,请大侠们帮忙。
回复
分享到:

使用道具 举报

 楼主| 发表于 2008-6-29 10:45 | 显示全部楼层
可能我的问题没有说清楚。我后面又在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.-')
发表于 2008-6-29 19:59 | 显示全部楼层

回复 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-')

评分

1

查看全部评分

 楼主| 发表于 2008-6-29 22:03 | 显示全部楼层

回复 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=[0.970:0.002:1,1.988:0.002:2];
spline_g_value=[nu_g(486:501),zeros(1,7)];
spline_f_value=[nu_f(486:501),zeros(1,7)];
spline_u_value=[nu_u(486:501),zeros(1,7)];

nuspline_t=0:0.002:2;
pp=spline(spline_t_value,spline_g_value);
nuspline_g=[nu_g,ppval(pp,1.002:0.002:2)];
pp=spline(spline_t_value,spline_f_value);
nuspline_f=[nu_f,ppval(pp,1.002:0.002:2)];
pp=spline(spline_t_value,spline_u_value);
nuspline_u=[nu_u,ppval(pp,1.002:0.002:2)];
normg=sqrt(intergral([0,tvalue],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),'-');
发表于 2008-6-30 17:26 | 显示全部楼层
不是很明白周期化延拓过程,这样做不是把原来的函数都变了么?

既然变了,那还有什么理由认为三者之间仍然满足原来的关系?



时域和频域的对应很简单的,

时域的采样点同频域的采样点个数相同,

时域两相邻采样点差的倒数为整个频域宽度,

同理,频域两相邻采样点差的倒数为整个时域宽度。



常用的代码如下:

ts=t(2)-t(1);

fs=1/ts;

f=(-N/2:N/2-1)/N*fs;
 楼主| 发表于 2008-6-30 18:07 | 显示全部楼层

回复 5楼 的帖子

恩,应该是你所说的这样。我也很迷惑。
但是文章上介绍的只是在t接近1时数据会受到影响。
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

QQ|小黑屋|Archiver|手机版|联系我们|声振论坛

GMT+8, 2024-9-23 07:24 , Processed in 0.057318 second(s), 18 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

快速回复 返回顶部 返回列表