微分几何,大哥们帮帮忙

(一)傅里叶反变换重建

傅里叶反变换重建法是一种重要的变换重建方法。它主要包括以下三个步骤：①建立数学模型，其中已和未知量都是连续实数的函数；②利用反变换公式解未知量；③调节反变换公式以适应离散、有噪声应用的需求。注意上面第二步中在解未知量时理论上可以有多个等价的公式来解。而第三步中，由于离散化时可采用不同的近似，所以理论上等价的方法对实际数据应用的结果会不同。在具体应用中，所测量到的数据对应于许多个离散点（s，θ）上的估计值g（s，θ），而所重建的图像也是一个离散的数组。

设在s和θ上都均匀采样的情况。考虑在N个相差△θ的角度上的投影，在森蠢每个角上用M个间距为△s的射线测量，定义整数M+和M-为：�

（6.181）

为了保证一系列射线｛（m△s，n△θ）：M�-≤m≤M+，1≤n≤N}覆盖单位圆，需要选△θ=π/N和△s=1／M+。此时g（m△s，n△θ）为平行投影的此毁陪射线数据。设图像区被一个直角网格覆盖，其中K+和K-用类似于式（6.181）的方法定义（K为X方向上的点数），L+和L-也用类似于式（6.181）的方法定义（L为Y方向上的点数）。根据这些定义，一个重建算法就是要通过M×N个测量值g（m△s，n△θ）估计出在K×L个采样点的f（K△x，l△y）。

1.傅里叶变换投影定理变换方法的基础是傅余侍里叶变换投影定理。设G（R，θ）是g（s，θ）对应第一个变量s的一维傅里叶变换，即

（6.182）

F（X，Y）是f（x，y）的二维傅里叶变换：

（6.183）

那么可以证明如下投影定理

G（R，θ）=F（Rcosθ，Rsinθ）（6.184）

即f（x，y）以θ角进行投影的傅里叶变换等于f（x，y）的傅里叶变换在傅里叶空间（R，θ）处的值。换句话说，f（x，y）在与X轴成θ角的直线上投影的傅里叶变换是f(x，y）的傅里叶变换在朝向角θ上的一个截面。

2.傅里叶反变换重建公式根据傅里叶变换投影定理很容易就可以得到傅里叶反变换重建的公式。对式（6.184）两边在直角坐标系中取傅里叶反变换

（6.185）

注意到G（·）是g（·）的傅里叶变换，所以上式是给出g（s，θ），计算f（x，y）的一个重建公式。这也是第一个在投影重建中得到应用的技术。实用中需要加一个窗，以把积分区限制在傅里叶XY平面上的一个有限区域W，这样得到f（x，y）的一个带限逼近fw（x，y）。

（6.186）

现在考虑计算G（·）。实际中G（·）只在一系列θ=θn角取值（θn代表nΔθ）。G（R，θn）可用在一系列采样点（mΔs，θn）对g（·）的求和得到。

如果令R=kΔR（k为整数，ΔR为采样），取ΔR=1／（MΔs），则

G∑（ΔR，θ△n）=Δs（mΔs，θn）exp〔-j2π（mΔs）〕（6.187）

根据G∑（kΔR，θn）可对任意的（X，Y）插值出G（X2＋Y2）1/2，arctgY\X。然后由式（6.186）可知

（6.188）

于是fw（x，y）可由下式确定

（6.189）

进一步，如果令△x=1/（U△X），△y=1／（V△Y），则

（6.190）

图6-134右半部给出在（傅里叶）XY平面上的（k△R，θn）位置，用“·”表示。在这些采样位置上可用式（6.187）计算G。这些点在平面上形成一个极坐标模式，各点在极坐标系中是等距分布的。左半部中“+”点表示位置（u△X，υ△Y），需要知道这些点的傅里叶变换以根据式（6.190）计算fw，（u△X，υ△Y）。这些点在平面上形成一个直角坐标网格。可根据在极坐标系中的已知值进行插值以得到这些网格点的估计值。上式只给出有限个（U个k，V个l）估计值fw（k△x，l△y）。对这些参数的双重求和可用FFT快速计算。图6-134（a）右半部的极坐标模式对傅里叶平面的采样不是很有效，一种可能的改进办法是采用图6-134（b）的模式。�

图6-134傅里叶空间的直角和极坐标网络

基于傅里叶反变换的重建技术主要有三个步骤：①对以角θn（n=1，2，…，N）的投影进行一维傅里叶变换；②在傅里叶空间从极坐标向直角坐标插值；③进行二维傅里叶反变换以得到重建图像，由于第三步需要用到二维变换，所以不能根据所获得的部分投影数据重建图像，而必须在获得全部投影数据后才能再重建图像。

在图象处理的广泛应用领域中，傅立叶变换起着非常重要的作用，具体表现在包括图象分析、图象增强及图象压缩等方面。

fftshift的作用正是让正半轴部分和负半轴部分的图像分别关于各自的中心对称。因为直接用fft得出的数据与频率不是对应的，fftshift可以纠正过来。

假设f（x,y）是一个离散空间中的二维函数，则该函数的二维傅立叶变换的定义如下：

p=0,1…M-1 q=0,1…N-1 (1)

或 p=0,1…M-1 q=0,1…N-1 （2）

离散傅立叶反变换的定义如下：

m=0,1…M-1 n=0,1…N-1（3）

F（p,q）称为f（m,n）的离散傅立叶变换系数。这个式子表明，函数f（m,n）可以用无数个不同频率的复指数信号和表示，而在频率（w1，w2）处的复指数信号的幅度和相位是F（w1，w2）。

2、MATLAB提供的快速傅立叶变换函数

（1）fft2

fft2函数用于计算二维快速傅立蠢备叶变换，其语法格式为：

B = fft2(I)

B = fft2(I)返回图象I的二维fft变换矩阵，输入图象I和输出图象B大小相同。

例如，计算图象的二维傅立叶变换，并显示其幅值的结果，如图所示，其命令格式如下

load imdemos saturn2

imshow(saturn2)

B = fftshift(fft2(saturn2))

imshow(log(abs(B)),[],'notruesize')

（2）fftshift

MATLAB提供的fftshift函数用于将变换后的图象频谱中心从矩阵的原点移到矩阵的中心，其语法格式为：

B = fftshift(I)

对于矩阵I，B = fftshift(I)将I的一、三象限和二、四象限进行互换。

（2）ifft2

ifft2函数用于计算图象的二维傅立叶反变换，其语法格睁缺式为：

B = ifft2(I)

B = ifft2(A)返回图象I的二维傅立叶反变换矩阵，输入图象I和输出图象B大小相同。其语法格式含义与fft2函数的语法格式相同，可以参考fft2函数的说明。

如果信号是二维的，用上面的函数即可！直接调用。

如果信号是一维的，给你下面的例子带早毁，你应该能明白！

clear

fs=100N=128 %采样频率和数据点数

n=0:N-1t=n/fs %时间序列

x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t)%信号

y=fft(x,N) %对信号进行快速Fourier变换，逆变换函数为ifft

mag=abs(y)%求得Fourier变换后的振幅

f=n*fs/N %频率序列

subplot(2,2,1),plot(f,mag) %绘出随频率变化的振幅

xlabel('频率/Hz')

ylabel('振幅')title('N=128')grid on

subplot(2,2,2),plot(f(1:N/2),mag(1:N/2))%绘出Nyquist频率之前随频率变化的振幅

xlabel('频率/Hz')

ylabel('振幅')title('N=128')grid on

%对信号采样数据为1024点的处理

fs=100N=1024n=0:N-1t=n/fs

x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t)%信号

y=fft(x,N) %对信号进行快速Fourier变换

mag=abs(y) %求取Fourier变换的振幅

f=n*fs/N

subplot(2,2,3),plot(f,mag)%绘出随频率变化的振幅

xlabel('频率/Hz')

ylabel('振幅')title('N=1024')grid on

subplot(2,2,4)

plot(f(1:N/2),mag(1:N/2))%绘出Nyquist频率之前随频率变化的振幅

xlabel('频率/Hz')

ylabel('振幅')title('N=1024')grid on

---