功率谱密度单位是什么？和功率有关系吗？
功率谱密度谱是一种概率统计方法，是对随机变量均方值的量度。一般用于随机振动分析，连续瞬态响应只能通过概率分布函数进行描述，即出现某水平响应所对应的概率。功率谱密度是结构在随机动态载荷激励下响应的统计结果，是一条功率谱密度值—频率值的关系曲线，其中功率谱密度可以是位移功率谱密度、速度功率谱密度、加速度功率谱密度、力功率谱密度等形式。数学上，功率谱密度值—频率值的关系曲线下的面积就是方差，即响应标准偏差的平方值。
谱是个很不严格的东西，常常指信号的Fourier变换，是一个时间平均(time average)概念功率谱的概念是针对功率有限信号的(能量有限信号可用能量谱分析)，所表现的是单位频带内信号功率随频率的变换情况。保留频谱的幅度信息，但是丢掉了相位信息，所以频谱不同的信号其功率谱是可能相同的。有两个重要区别： 1。功率谱是随机过程的统计平均概念，平稳随机过程的功率谱是一个确定函数;而频谱是随机过程样本的Fourier变换，对于一个随机过程而言，频谱也是一个“随机过程”。(随机的频域序列) 2。功率概念和幅度概念的差别。此外，只能对宽平稳的各态历经的二阶矩过程谈功率谱，其存在性取决于二阶局是否存在并且二阶矩的Fourier变换收敛;而频谱的存在性仅仅取决于该随机过程的该样本的Fourier变换是否收敛。热心网友回答提问者对于答案的评价：谢谢解答。
频谱分析(也称频率分析)，是对动态信号在频率域内进行分析，分析的
结果是以频率为坐标的各种物理量的谱线和曲线，可得到各种幅值以频率为变
量的频谱函数F(ω)。频谱分析中可求得幅值谱、相位谱、功率谱和各种谱密
度等等。频谱分析过程较为复杂，它是以傅里叶级数和傅里叶积分为基础的。
功率谱是个什么概念?它有单位吗?
随机信号是时域无限信号，不具备可积分条件，因此不能直接进行傅氏变换。一般用具有统计特性的功率谱来作为谱分析的依据。功率谱与自相关函数是一个傅氏变换对。功率谱具有单位频率的平均功率量纲。所以标准叫法是功率谱密度。通过功率谱密度函数，可以看出随机信号的能量随着频率的分布情况。像白噪声就是平行于w轴，在w轴上方的一条直线。
功率谱密度，从名字分解来看就是说，观察对象是功率，观察域是谱域，通常指频域，密度，就是指观察对象在观察域上的分布情况。一般我们讲的功率谱密度都是针对平稳随机过程的，由于平稳随机过程的样本函数一般不是绝对可积的，因此不能直接对它进行傅立叶分析。可以有三种办法来重新定义谱密度，来克服上述困难。
一是用相关函数的傅立叶变换来定义谱密度;二是用随机过程的有限时间傅立叶变换来定义谱密度;三是用平稳随机过程的谱分解来定义谱密度。三种定义方式对应于不同的用处，首先第一种方式前提是平稳随机过程不包含周期分量并且均值为零，这样才能保证相关函数在时差趋向于无穷时衰减，所以lonelystar说的不全对，光靠相关函数解决不了许多问题，要求太严格了;对于第二种方式，虽然一个平稳随机过程在无限时间上不能进行傅立叶变换，但是对于有限区间，傅立叶变换总是存在的，可以先架构有限时间区间上的变换，在对时间区间取极限，这个定义方式就是当前快速傅立叶变换(FFT)估计谱密度的依据;第三种方式是根据维纳的广义谐和分析理论：Generalized harmonic analysis, Acta Math, 55(1930),117-258,利用傅立叶-斯蒂吉斯积分，对均方连续的零均值平稳随机过程进行重构，在依靠正交性来建立的。
另外，对于非平稳随机过程，也有三种谱密度建立方法，由于字数限制，功率谱密度的单位是G的平方/频率。就是就是函数幅值的均方根值与频率之比。是对随机振动进行分析的重要参数。
功率谱密度的国际单位是什么?
如果是加速度功率谱密度，加速度的单位是m/s^2,
那么，加速度功率谱密度的单位就是(m/s^2)^2/Hz,
而Hz的单位是1/s,经过换算得到加速度功率谱密度的单位是m^2/s^3.
同理，如果是位移功率谱密度，它的单位就是m^2*s,
如果是弯矩功率谱密度，单位就是(N*m)^2*s
位移功率谱——m^2*s
速度功率谱——m^2/s
加速度功率谱——m^2/s^3
信号傅立叶变换的幅度图和频谱图matlab示例
fs=100;
x=-2:1/fs:2;
y=sin(3*pi*x);
z=rectpuls(x);
figure;plot(x,y,x,z,':r');
my=abs(fft(y));
mz=abs(fft(z));
my=my/max(my); %归一化
mz=mz/max(mz); %归一化
f=(0:1/length(x):1)*fs;
figure;plot(f(1:fs/2),my(1:fs/2),f(1:fs/2),mz(1:fs/2),':r');
my=20*log10(my+eps);mz=20*log10(mz+eps);
figure;plot(f(1:fs/2),my(1:fs/2),f(1:fs/2),mz(1:fs/2),':r'
时域信号--->相关函数--(FFT变换)-->功率谱--(除以频率分辨率)-->功率谱密度，这叫做间接求法，可以抑制白噪声，或者通俗的说不规律信号，分析的点数越多，规律信号的信噪比越好。
时域信号--(FFT变换)-->幅度谱--(平方)-->功率谱，这叫直接求法，最好不要用，除非你就想分析噪声有多大
转自：http://wenda.tianya.cn/question/4c69f62914e8403d
MATLAB中FFT的使用方法
说明：以下资源来源于《数字信号处理的MATLAB实现》万永革主编一.调用方法
X=FFT(x);
X=FFT(x，N);
x=IFFT(X);
x=IFFT(X,N)
用MATLAB进行谱分析时注意：
(1)函数FFT返回值的数据结构具有对称性。
例：
N=8;
n=0:N-1;
xn=;
Xk=fft(xn)
→
Xk =
39.0000 -10.7782 + 6.2929i 0 - 5.0000i 4.7782 - 7.7071i 5.0000 4.7782 + 7.7071i 0 + 5.0000i -10.7782 - 6.2929i
Xk与xn的维数相同，共有8个元素。Xk的第一个数对应于直流分量，即频率值为0。
(2)做FFT分析时，幅值大小与FFT选择的点数有关，但不影响分析结果。在IFFT时已经做了处理。要得到真实的振幅值的大小，只要将得到的变换后结果乘以2除以N即可。
二.FFT应用举例
例1：x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t)。采样频率fs=100Hz，分别绘制N=128、1024点幅频图。
clf;
fs=100;N=128; %采样频率和数据点数
n=0:N-1;t=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,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=100;N=1024;n=0:N-1;t=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;
运行结果：
fs=100Hz，Nyquist频率为fs/2=50Hz。整个频谱图是以Nyquist频率为对称轴的。并且可以明显识别出信号中含有两种频率成分：15Hz和40Hz。由此可以知道FFT变换数据的对称性。因此用FFT对信号做谱分析，只需考察0~Nyquist频率范围内的福频特性。若没有给出采样频率和采样间隔，则分析通常对归一化频率0~1进行。另外，振幅的大小与所用采样点数有关，采用128点和1024点的相同频率的振幅是有不同的表现值，但在同一幅图中，40Hz与15Hz振动幅值之比均为4：1，与真实振幅0.5：2是一致的。为了与真实振幅对应，需要将变换后结果乘以2除以N。
例2：x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t),fs=100Hz,绘制：
(1)数据个数N=32，FFT所用的采样点数NFFT=32;
(2)N=32，NFFT=128;
(3)N=136，NFFT=128;
(4)N=136，NFFT=512。
clf;fs=100; %采样频率
Ndata=32; %数据长度
N=32; %FFT的数据长度
n=0:Ndata-1;t=n/fs; %数据对应的时间序列
x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t); %时间域信号
y=fft(x,N); %信号的Fourier变换
mag=abs(y); %求取振幅
f=(0:N-1)*fs/N; %真实频率
subplot(2,2,1),plot(f(1:N/2),mag(1:N/2)*2/N); %绘出Nyquist频率之前的振幅
xlabel('频率/Hz');ylabel('振幅');
title('Ndata=32 Nfft=32');grid on;
Ndata=32; %数据个数
N=128; %FFT采用的数据长度
n=0:Ndata-1;t=n/fs; %时间序列
x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t);
y=fft(x,N);
mag=abs(y);
f=(0:N-1)*fs/N; %真实频率
subplot(2,2,2),plot(f(1:N/2),mag(1:N/2)*2/N); %绘出Nyquist频率之前的振幅
xlabel('频率/Hz');ylabel('振幅');
title('Ndata=32 Nfft=128');grid on;
Ndata=136; %数据个数
N=128; %FFT采用的数据个数
n=0:Ndata-1;t=n/fs; %时间序列
x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t);
y=fft(x,N);
mag=abs(y);
f=(0:N-1)*fs/N; %真实频率
subplot(2,2,3),plot(f(1:N/2),mag(1:N/2)*2/N); %绘出Nyquist频率之前的振幅
xlabel('频率/Hz');ylabel('振幅');
title('Ndata=136 Nfft=128');grid on;
Ndata=136; %数据个数
N=512; %FFT所用的数据个数
n=0:Ndata-1;t=n/fs; %时间序列
x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t);
y=fft(x,N);
mag=abs(y);
f=(0:N-1)*fs/N; %真实频率
subplot(2,2,4),plot(f(1:N/2),mag(1:N/2)*2/N); %绘出Nyquist频率之前的振幅
xlabel('频率/Hz');ylabel('振幅');
title('Ndata=136 Nfft=512');grid on;
结论：
(1)当数据个数和FFT采用的数据个数均为32时，频率分辨率较低，但没有由于添零而导致的其他频率成分。
(2)由于在时间域内信号加零，致使振幅谱中出现很多其他成分，这是加零造成的。其振幅由于加了多个零而明显减小。
(3)FFT程序将数据截断，这时分辨率较高。
(4)也是在数据的末尾补零，但由于含有信号的数据个数足够多，FFT振幅谱也基本不受影响。
对信号进行频谱分析时，数据样本应有足够的长度，一般FFT程序中所用数据点数与原含有信号数据点数相同，这样的频谱图具有较高的质量，可减小因补零或截断而产生的影响。
例3：x=cos(2*pi*0.24*n)+cos(2*pi*0.26*n)
(1)数据点过少，几乎无法看出有关信号频谱的详细信息;
(2)中间的图是将x(n)补90个零，幅度频谱的数据相当密，称为高密度频谱图。但从图中很难看出信号的频谱成分。
(3)信号的有效数据很长，可以清楚地看出信号的频率成分，一个是0.24Hz，一个是0.26Hz，称为高分辨率频谱。
可见，采样数据过少，运用FFT变换不能分辨出其中的频率成分。添加零后可增加频谱中的数据个数，谱的密度增高了，但仍不能分辨其中的频率成分，即谱的分辨率没有提高。只有数据点数足够多时才能分辨其中的频率成分。
转自：http://blog.163.com/fei_lai_feng/blog/static/9289962200971751114547/
理解快速傅里叶变换（FFT）算法
快速傅里叶变换(Fast Fourier Transform)是信号处理与数据分析领域里最重要的算法之一。没有正规计算机科学课程背景的我，使用这个算法多年，但这周我却突然想起自己从没思考过为什么FFT能如此快速地计算离散傅里叶变换。我打开一本老旧的算法书，欣赏了JW Cooley 和 John Tukey 在1965年的文章中，以看似简单的计算技巧来讲解这个东西。本文的目标是，深入Cooley-Tukey FFT 算法，解释作为其根源的“对称性”，并以一些直观的python代码将其理论转变为实际。我希望这次研究能使数据科学家(例如我)，对这个算法的背景原理有更全面的认识。
FFT(快速傅里叶变换)本身就是离散傅里叶变换(Discrete Fourier Transform)的快速算法，使算法复杂度由原本的O(N^2) 变为 O(NlogN)，离散傅里叶变换DFT，如同更为人熟悉的连续傅里叶变换，有如下的正、逆定义形式：
xn 到 Xk 的转化就是空域到频域的转换，这个转换有助于研究信号的功率谱，和使某些问题的计算更有效率。事实上，你还可以查看一下我们即将推出的天文学和统计学的图书的第十章(这里有一些图示和python代码)。作为一个例子，你可以查看下我的文章《用python求解薛定谔方程》，是如何利用FFT将原本复杂的微分方程简化。
正因为FFT在那么多领域里如此有用，python提供了很多标准工具和封装来计算它。NumPy 和 SciPy 都有经过充分测试的封装好的FFT库，分别位于子模块 numpy.fft 和 scipy.fftpack 。我所知的最快的FFT是在 FFTW包中 ，而你也可以在python的pyFFTW 包中使用它。
虽然说了这么远，但还是暂时先将这些库放一边，考虑一下怎样使用原始的python从头开始计算FFT。
计算离散傅里叶变换
简单起见，我们只关心正变换，因为逆变换也只是以很相似的方式就能做到。看一下上面的DFT表达式，它只是一个直观的线性运算：向量x的矩阵乘法，
矩阵M可以表示为
这么想的话，我们可以简单地利用矩阵乘法计算DFT：
import numpy as np
def DFT_slow(x):
"""Compute the discrete Fourier Transform of the 1D array x"""
x = np.asarray(x, dtype=float)
N = x.shape
n = np.arange(N)
k = n.reshape((N, 1))
M = np.exp(-2j * np.pi * k * n / N)
return np.dot(M, x)
对比numpy的内置FFT函数，我们来对结果进行仔细检查，
x = np.random.random(1024)
np.allclose(DFT_slow(x), np.fft.fft(x))
输出：
True
现在为了验证我们的算法有多慢，对比下两者的执行时间
%timeit DFT_slow(x)
%timeit np.fft.fft(x)
输出：
10 loops, best of 3: 75.4 ms per loop
10000 loops, best of 3: 25.5 μs per loop
使用这种简化的实现方法，正如所料，我们慢了一千多倍。但问题不是这么简单。对于长度为N的输入矢量，FFT是O(N logN)级的，而我们的慢算法是O(N^2)级的。这就意味着，FFT用50毫秒能干完的活，对于我们的慢算法来说，要差不多20小时! 那么FFT是怎么提速完事的呢?答案就在于他利用了对称性。
离散傅里叶变换中的对称性
算法设计者所掌握的最重要手段之一，就是利用问题的对称性。如果你能清晰地展示问题的某一部分与另一部分相关，那么你就只需计算子结果一次，从而节省了计算成本。
Cooley 和 Tukey 正是使用这种方法导出FFT。 首先我们来看下
的值。根据上面的表达式，推出：
对于所有的整数n，exp=1。
最后一行展示了DFT很好的对称性：
简单地拓展一下：
同理对于所有整数 i 。正如下面即将看到的，这个对称性能被利用于更快地计算DFT。
DFT 到 FFT：
利用对称性 Cooley 和 Tukey 证明了，DFT的计算可以分为两部分。从DFT的定义得：
我们将单个DFT分成了看起来相似两个更小的DFT。一个奇，一个偶。目前为止，我们还没有节省计算开销，每一部分都包含(N/2)∗N的计算量，总的来说，就是N^2 。
技巧就是对每一部分利用对称性。因为 k 的范围是 0≤k
但我们不是到这步就停下来，只要我们的小傅里叶变换是偶倍数，就可以再作分治，直到分解出来的子问题小到无法通过分治提高效率，接近极限时，这个递归是 O(n logn) 级的。
这个递归算法能在python里快速实现，当子问题被分解到合适大小时，再用回原本那种“慢方法”。
def FFT(x):
"""A recursive implementation of the 1D Cooley-Tukey FFT"""
x = np.asarray(x, dtype=float)
N = x.shape
if N % 2 > 0:
raise ValueError("size of x must be a power of 2")
elif N <= 32:
# this cutoff should be optimized
return DFT_slow(x)
else:
X_even = FFT(x[::2])
X_odd = FFT(x)
factor = np.exp(-2j * np.pi * np.arange(N) / N)
return np.concatenate( * X_odd,
X_even + factor * X_odd])
现在我们做个快速的检查，看结果是否正确：
x = np.random.random(1024)
np.allclose(FFT(x), np.fft.fft(x))
True
然后与“慢方法”的运行时间对比下：
%timeit DFT_slow(x)
%timeit FFT(x)
%timeit np.fft.fft(x)
10 loops, best of 3: 77.6 ms per loop
100 loops, best of 3: 4.07 ms per loop
10000 loops, best of 3: 24.7 μs per loop
现在的算法比之前的快了一个数量级。而且，我们的递归算法渐近于 O(n logn) 。我们实现了FFT 。
需要注意的是，我们还没做到numpy的内置FFT算法，这是意料之中的。numpy 的 fft 背后的FFTPACK算法 是以 Fortran 实现的，经过了多年的调优。此外，我们的NumPy的解决方案，同时涉及的Python堆栈递归和许多临时数组的分配，这显著地增加了计算时间。
还想加快速度的话，一个好的方法是使用Python/ NumPy的工作时，尽可能将重复计算向量化。我们是可以做到的，在计算过程中消除递归，使我们的python FFT更有效率。
向量化的NumPy
注意上面的递归FFT实现，在最底层的递归，我们做了N/32次的矩阵向量乘积。我们的算法会得益于将这些矩阵向量乘积化为一次性计算的矩阵-矩阵乘积。在每一层的递归，重复的计算也可以被向量化。因为NumPy很擅长这类操作，我们可以利用这一点来实现向量化的FFT。
def FFT_vectorized(x):
"""A vectorized, non-recursive version of the Cooley-Tukey FFT"""
x = np.asarray(x, dtype=float)
N = x.shape
if np.log2(N) % 1 > 0:
raise ValueError("size of x must be a power of 2")
# N_min here is equivalent to the stopping condition above,
# and should be a power of 2
N_min = min(N, 32)
# Perform an O DFT on all length-N_min sub-problems at once
n = np.arange(N_min)
k = n[:, None]
M = np.exp(-2j * np.pi * n * k / N_min)
X = np.dot(M, x.reshape((N_min, -1)))
# build-up each level of the recursive calculation all at once
while X.shape < N:
X_even = X[:, :X.shape / 2]
X_odd = X[:, X.shape / 2:]
factor = np.exp(-1j * np.pi * np.arange(X.shape)
/ X.shape)[:, None]
X = np.vstack([X_even + factor * X_odd,
X_even - factor * X_odd])
return X.ravel()
x = np.random.random(1024)
np.allclose(FFT_vectorized(x), np.fft.fft(x))
True
因为我们的算法效率更大幅地提升了，所以来做个更大的测试(不包括DFT_slow)
x = np.random.random(1024 * 16)
%timeit FFT(x)
%timeit FFT_vectorized(x)
%timeit np.fft.fft(x)
10 loops, best of 3: 72.8 ms per loop
100 loops, best of 3: 4.11 ms per loop
1000 loops, best of 3: 505 μs per loop
我们的实现又提升了一个级别。这里我们是以 FFTPACK中大约10以内的因数基准，用了仅仅几十行 Python + NumPy代码。虽然没有相应的计算来证明， Python版本是远优于 FFTPACK源，这个你可以从这里浏览到。
那么 FFTPACK是怎么获得这个最后一点的加速的呢?也许它只是一个详细的记录簿， FFTPACK花了大量时间来保证任何的子计算能够被复用。我们这里的numpy版本涉及到额外的内存的分配和复制，对于如Fortran的一些低级语言就能够很容易的控制和最小化内存的使用。并且Cooley-Tukey算法还能够使其分成超过两部分(正如我们这里用到的Cooley-Tukey FFT基2算法)，而且，其它更为先进的FFT算法或许也可以能够得到应用，包括基于卷积的从根本上不同的方法(例如Bluestein的算法和Rader的算法)。结合以上的思路延伸和方法，就可使阵列大小即使不满足2的幂，FFT也能快速执行。
我希望他们提供大量的基于FFT数据分析是怎样进行的背景，尽管纯Python函数在实际中并不适用。作为数据科学家，我们可以暂且引进能够包含基本应用工具的黑盒子，是由我们有更多算法思想的同事参入所构建，但是我坚定的相信：当我们对所应用的数据的底层算法有更深的理解，我们也将会成为更优秀的从业者。
转自：http://blog.jobbole.com/58246/
FFT。。。。。。。
快速傅氏变换，是离散傅氏变换的快速算法，它是根据离散傅氏变换的奇、偶、虚、实等特性，对离散傅立叶变换的算法进行改进获得的。它对傅氏变换的理论并没有新的发现，但是对于在计算机系统或者说数字系统中应用离散傅立叶变换，可以说是进了一大步。 设x(n)为N项的复数序列，由DFT变换，任一X(m)的计算都需要N次复数乘法和N-1次复数加法，而一次复数乘法等于四次实数乘法和两次实数加法，一次复数加法等于两次实数加法，即使把一次复数乘法和一次复数加法定义成一次“运算”(四次实数乘法和四次实数加法)，那么求出N项复数序列的X(m), 即N点DFT变换大约就需要N2次"运算"。当N=1024点甚至更多的时候，需要N2=1048576次运算.在FFT中，利用WN的周期性和对称性，把一个N项序列(设N=2k,k为正整数)，分为两个N/2项的子序列，每个N/2点DFT变换需要(N/2)2次运算，再用N次运算把两个N/2点的DFT 变换组合成一个N点的DFT变换。这样变换以后，总的运算次数就变成N+2(N/2)2=N+N2/2。继续上面的例子，N=1024时，总的运算次数就变成了525312次，节省了大约50%的运算量。而如果我们将这种“一分为二”的思想不断进行下去，直到分成两两一组的DFT运算单元，那么N点的 DFT变换就只需要Nlog2N次的运算，N在1024点时，运算量仅有10240次，是先前的直接算法的1%，点数越多，运算量的节约就越大，这就是 FFT的优越性. 傅里叶变换(Transformée de Fourier)是一种积分变换。因其基本思想首先由法国学者傅里叶系统地提出，所以以其名字来命名以示纪念。 应用 傅里叶变换在物理学、数论、组合数学、信号处理、概率论、统计学、密码学、声学、光学、海洋学、结构动力学等领域都有着广泛的应用(例如在信号处理中，傅里叶变换的典型用途是将信号分解成幅值分量和频率分量)。 概要介绍 傅里叶变换能将满足一定条件的某个函数表示成三角函数(正弦和/或余弦函数)或者它们的积分的线性组合。在不同的研究领域，傅里叶变换具有多种不同的变体形式，如连续傅里叶变换和离散傅里叶变换。最初傅里叶分析是作为热过程的解析分析的工具被提出的(参见：林家翘、西格尔著《自然科学中确定性问题的应用数学》，科学出版社，北京。原版书名为 C. C. Lin & L. A. Segel, Mathematics Applied to Deterministic Problems in the Natural Sciences, Macmillan Inc., New York, 1974)。 傅里叶变换属于谐波分析。 傅里叶变换的逆变换容易求出,而且形式与正变换非常类似; 正弦基函数是微分运算的本征函数,从而使得线性微分方程的求解可以转化为常系数的代数方程的求解.在线性时不变的物理系统内,频率是个不变的性质,从而系统对于复杂激励的响应可以通过组合其对不同频率正弦信号的响应来获取; 卷积定理指出:傅里叶变换可以化复杂的卷积运算为简单的乘积运算,从而提供了计算卷积的一种简单手段; 离散形式的傅里叶变换可以利用数字计算机快速的算出(其算法称为快速傅里叶变换算法(FFT)).
转自：http://blog.sina.com.cn/s/blog_60484b020100dcjd.html
[FFT] FFT初解 。。。。。
FFT初解谨以此文纪念过往的岁月
一.前言
首先申明俺不是一个算法工程师，俺是一个底层驱动工程师，有人会发问一个底层驱动工程师需要这个吗?但是我不幸的告诉你，确实是需要的，不过我们不要像算法工程师那样搞得很精通，但是还是需要去了解这是个什么东西。说实话，这个东西在大学时候学过，还好好的去理解了一样，不过到现在忘的差不多了，这愈发的让我明白一句话，好记性不如烂笔头，如果以前有好好记录的好习惯，那现在只要把以前的东西拿出来看看再印证一下就可以了。不过历史不可以如果。为了不让明天继续懊悔今天，在这里记录下本人学习的一些记录。
二.FFT的数学基础
FFT是什么，额，说白了就是一种时域和频域变换的手段。什么是时域什么是频域，额，你google吧，鄙视百度!
这里说明一下，相比较FFT而言，我们更加关心离散傅里叶变换(DTF),因为我们底层驱动工程师面对的是一个一个离散的数据。也许大家还是对FFT的变换的概念还不理解，那我们举个例子来说吧：
1 AD在一个1Hz正弦波周期采集了1024个数据(即你的采样频率为1024Hz)，从时域上看1s内采集了1024数据，那1024个数据做1024点的FFT变换。我们以数据的下标作为横轴，而数据的值作为纵轴，你会发现第一个点的值最大，我们可以从该值计算出我们被采样的频率为1Hz，值是这样算出来的1024(Hz)/1024*1 = 1Hz,如果最大的值点是第8个点，那被采样频率是8Hz。不过要满足乃奎斯定律，即采样频率大于被采样频率的两倍。
为什么会有东西出现呢，主要是时域的信号好采样，至于频域的信号难以采样，于是傅里叶这个天才就发明了这个东东，不过这可是现代数字信号处理的基础。
闲话不说了，还是来看DFT的数学基础(下面的内容不少拷贝wiki的)：
对于复数序列
，离散傅里叶变换公式为：
下面，我们用N次单位根WN来表示
。
WN的性质：
周期性，WN具有周期N。
对称性：
。
为了简单起见，我们下面设待变换序列长度n = 2r。 根据上面单位根的对称性，求级数
时，可以将求和区间分为两部分：
Fodd(k) 和 Feven(k)是两个分别关于序列
奇数号和偶数号序列N/2点变换。由此式只能计算出yk的前N/2个点，对于后N/2个点，注意 Fodd(k) 和 Feven(k) 都是周期为N/2的函数，由单位根的对称性，于是有以下变换公式：
。
这样，一个N点变换就分解成了两个N/2点变换。照这样可继续分解下去。这就是库利-图基快速傅里叶变换算法的基本原理。根据主定理不难分析出此时算法的时间复杂度为O(NlogN)
三.8点快速傅里叶变换
主要是蝶形算法
FFT算法图
__no_init float fft_r; //FFT输入实数序列，保存变换后的频域实部
__no_init float fft_i; //保存变换后的频域虚部
__no_init float harm; //各次谐波幅值
const Uint8 FFT_BIT={ //供完成倒位序排列使用的位定义
0x01,0x02,0x04,0x08,0x10,0x20,0x40,0x80};
void fft(void)
{
Uint16 i,j,k,L;
Uint16 b,c,p,n;
Uint8 x,xc;
float Gr,Gi,Wr,Wi; for(i=0;i
{
x =i;
xc=0;
for(j=0;j
{
if((x&0x01)!=0)
{
xc|=FFT_BIT[(FFT_M-1)-j];
}
x>>=1;
}
fft_i=fft_r;
}
for(i=0;i
{
fft_r=fft_i;
fft_i=0;
}
//以上是倒叙就是将数组中的数据进行倒叙，即将数组序号倒序，如FFT_BIT存储到FFT_BIT,为什么呢?5 = 0x110，进行二进制倒转变为0x011,变为3.
for(L=1;L<=FFT_M;L++)
{
b=1;
p=FFT_N>>1;
for(i=0;i
{
b<<=1; //b=2^(L-1)
p>>=1; //p=N/(2^L)
}
c=b<<1; //c=b*2
for(j=0;j
{
n=p*j;
for(k=j;k
{
Gr = fft_r;
Gi = fft_i;
Wr = Sinf*fft_r + Sinf*fft_i;
Wi = Sinf*fft_i - Sinf*fft_r;
fft_r = Gr + Wr;
fft_i = Gi + Wi;
fft_r = Gr - Wr;
fft_i = Gi - Wi;
}
}
}
}
//以上循环分为三个一个根据变换点的阶数循环，即上图中level1->level2->level3。第二个循环是group步进，level1为2 level 2为4。 第三个循环是需要两两相乘数据之间的length，如level1的length为1，level2的length为2，level3的length为4。
关于旋转因子，就是cos和sin的值，不过这两者之间是可以互相转换的。
也许到此大家对于快速傅里叶有了大概的理解，下面数C++的code，也是源于网上。这段code更好的描述的FFT的蝶形算法，值得参考。
void FFT(unsigned long & ulN, vector >& vecList)
{
//得到幂数
unsigned long ulPower = 0; //幂数
unsigned long ulN1 = ulN - 1;
while(ulN1 > 0)
{
ulPower++;
ulN1 /= 2;
}
//反序
bitset bsIndex; //二进制容器
unsigned long ulIndex; //反转后的序号
unsigned long ulK;
for(unsigned long p = 0; p < ulN; p++)
{
ulIndex = 0;
ulK = 1;
bsIndex = bitset(p);
for(unsigned long j = 0; j < ulPower; j++)
{
ulIndex += bsIndex.test(ulPower - j - 1) ? ulK : 0;
ulK *= 2;
}
if(ulIndex > p)
{
complex c = vecList;
vecList = vecList;
vecList = c;
}
}
//计算旋转因子
vector > vecW;
for(unsigned long i = 0; i < ulN / 2; i++)
{
vecW.push_back(complex(cos(2 * i * PI / ulN) , -1 * sin(2 * i * PI / ulN)));
}
for(unsigned long m = 0; m < ulN / 2; m++)
{
cout<< "\nvW[" << m << "]=" << vecW;
}
//计算FFT
unsigned long ulGroupLength = 1; //段的长度
unsigned long ulHalfLength = 0; //段长度的一半
unsigned long ulGroupCount = 0; //段的数量
complex cw; //WH(x)
complex c1; //G(x) + WH(x)
complex c2; //G(x) - WH(x)
for(unsigned long b = 0; b < ulPower; b++)
{
ulHalfLength = ulGroupLength;
ulGroupLength *= 2;
for(unsigned long j = 0; j < ulN; j += ulGroupLength)
{
for(unsigned long k = 0; k < ulHalfLength; k++)
{
cw = vecW * vecList;
c1 = vecList + cw;
c2 = vecList - cw;
vecList = c1;
vecList = c2;
}
}
}
}
四.总结
FFT数字信号处理的重中之重。
转自：http://blog.chinaunix.net/uid-24517893-id-336736.html
不要粘来的东西在论坛里胡搞。建议删除！
快速傅立叶变换的问题
快速傅氏变换，是离散傅氏变换的快速算法，它是根据离散傅氏变换的奇、偶、虚、实等特性，对离散傅立叶变换的算法进行改进获得的。它对傅氏变换的理论并没有新的发现，但是对于在计算机系统或者说数字系统中应用离散傅立叶变换，可以说是进了一大步。设x(n)为N项的复数序列，由DFT变换，任一X(m)的计算都需要N次复数乘法和N-1次复数加法，而一次复数乘法等于四次实数乘法和两次实数加法，一次复数加法等于两次实数加法，即使把一次复数乘法和一次复数加法定义成一次“运算”(四次实数乘法和四次实数加法)，那么求出N项复数序列的X(m), 即N点DFT变换大约就需要N2次运算。当N=1024点甚至更多的时候，需要N2=1048576次运算，在FFT中，利用WN的周期性和对称性，把一个N项序列(设N=2k,k为正整数)，分为两个N/2项的子序列，每个N/2点DFT变换需要(N/2)2次运算，再用N次运算把两个N/2点的DFT 变换组合成一个N点的DFT变换。这样变换以后，总的运算次数就变成N+2(N/2)2=N+N2/2。继续上面的例子，N=1024时，总的运算次数就变成了525312次，节省了大约50%的运算量。而如果我们将这种“一分为二”的思想不断进行下去，直到分成两两一组的DFT运算单元，那么N点的 DFT变换就只需要Nlog2N次的运算，N在1024点时，运算量仅有10240次，是先前的直接算法的1%，点数越多，运算量的节约就越大，这就是 FFT的优越性.
傅里叶变换(Transformée de Fourier)是一种积分变换。因其基本思想首先由法国学者傅里叶系统地提出，所以以其名字来命名以示纪念。
应用
傅里叶变换在物理学、数论、组合数学、信号处理、概率论、统计学、密码学、声学、光学、海洋学、结构动力学等领域都有着广泛的应用(例如在信号处理中，傅里叶变换的典型用途是将信号分解成幅值分量和频率分量)。
概要介绍
傅里叶变换能将满足一定条件的某个函数表示成三角函数(正弦和/或余弦函数)或者它们的积分的线性组合。在不同的研究领域，傅里叶变换具有多种不同的变体形式，如连续傅里叶变换和离散傅里叶变换。最初傅里叶分析是作为热过程的解析分析的工具被提出的(参见：林家翘、西格尔著《自然科学中确定性问题的应用数学》，科学出版社，北京。原版书名为 C. C. Lin & L. A. Segel, Mathematics Applied to Deterministic Problems in the Natural Sciences, Macmillan Inc., New York, 1974)。
傅里叶变换属于谐波分析。
傅里叶变换的逆变换容易求出,而且形式与正变换非常类似;
正弦基函数是微分运算的本征函数,从而使得线性微分方程的求解可以转化为常系数的代数方程的求解.在线性时不变的物理系统内,频率是个不变的性质,从而系统对于复杂激励的响应可以通过组合其对不同频率正弦信号的响应来获取;
卷积定理指出:傅里叶变换可以化复杂的卷积运算为简单的乘积运算,从而提供了计算卷积的一种简单手段;
离散形式的傅里叶变换可以利用数字计算机快速的算出(其算法称为快速傅里叶变换算法(FFT)).
转自：http://blog.sina.com.cn/s/blog_4b2f8f9601008zqh.html
傅立叶变换在图像处理中的作用
本帖最后由 wdhd 于 2016-3-21 11:03 编辑从现代数学的眼光来看，傅里叶变换是一种特殊的积分变换。它能将满足一定条件的某个函数表示成正弦基函数的线性组合或者积分。在不同的研究领域，傅里叶变换具有多种不同的变体形式，如连续傅里叶变换和离散傅里叶变换。
傅立叶变换属于调和分析的内容。"分析"二字，可以解释为深入的研究。从字面上来看，"分析"二字，实际就是"条分缕析"而已。它通过对函数的"条分缕析"来达到对复杂函数的深入理解和研究。从哲学上看，"分析主义"和"还原主义"，就是要通过对事物内部适当的分析达到增进对其本质理解的目的。比如近代原子论试图把世界上所有物质的本源分析为原子，而原子不过数百种而已，相对物质世界的无限丰富，这种分析和分类无疑为认识事物的各种性质提供了很好的手段。
在数学领域，也是这样，尽管最初傅立叶分析是作为热过程的解析分析的工具，但是其思想方法仍然具有典型的还原论和分析主义的特征。"任意"的函数通过一定的分解，都能够表示为正弦函数的线性组合的形式，而正弦函数在物理上是被充分研究而相对简单的函数类，这一想法跟化学上的原子论想法何其相似!奇妙的是,现代数学发现傅立叶变换具有非常好的性质,使得它如此的好用和有用,让人不得不感叹造物的神奇:
1. 傅立叶变换是线性算子,若赋予适当的范数,它还是酉算子;
2. 傅立叶变换的逆变换容易求出,而且形式与正变换非常类似;
3. 正弦基函数是微分运算的本征函数,从而使得线性微分方程的求解可以转化为常系数的代数方程的求解.在线性时不变的物理系统内,频率是个不变的性质,从而系统对于复杂激励的响应可以通过组合其对不同频率正弦信号的响应来获取;
4. 著名的卷积定理指出:傅立叶变换可以化复杂的卷积运算为简单的乘积运算,从而提供了计算卷积的一种简单手段;
5. 离散形式的傅立叶变换可以利用数字计算机快速的算出(其算法称为快速傅立叶变换算法(FFT)).
正是由于上述的良好性质,傅里叶变换在物理学、数论、组合数学、信号处理、概率、统计、密码学、声学、光学等领域都有着广泛的应用。
傅立叶变换在图像处理中有非常非常的作用
傅立叶变换在图像处理中有非常非常的作用。因为不仅傅立叶分析涉及图像处理的很多方面，傅立叶的改进算法，
比如离散余弦变换，gabor与小波在图像处理中也有重要的分量。
印象中，傅立叶变换在图像处理以下几个话题都有重要作用：
1.图像增强与图像去噪
绝大部分噪音都是图像的高频分量，通过低通滤波器来滤除高频——噪声; 边缘也是图像的高频分量，可以通过添加高频分量来增强原始图像的边缘;
2.图像分割之边缘检测
提取图像高频分量
3.图像特征提取：
形状特征：傅里叶描述子
纹理特征：直接通过傅里叶系数来计算纹理特征
其他特征：将提取的特征值进行傅里叶变换来使特征具有平移、伸缩、旋转不变性
4.图像压缩
可以直接通过傅里叶系数来压缩数据;常用的离散余弦变换是傅立叶变换的实变换;
傅立叶变换
傅里叶变换是将时域信号分解为不同频率的正弦信号或余弦函数叠加之和。连续情况下要求原始信号在一个周期内满足绝对可积条件。离散情况下，傅里叶变换一定存在。冈萨雷斯版<图像处理>里面的解释非常形象：一个恰当的比喻是将傅里叶变换比作一个玻璃棱镜。棱镜是可以将光分解为不同颜色的物理仪器，每个成分的颜色由波长(或频率)来决定。傅里叶变换可以看作是数学上的棱镜，将函数基于频率分解为不同的成分。当我们考虑光时,讨论它的光谱或频率谱。同样,傅立叶变换使我们能通过频率成分来分析一个函数。
傅立叶变换有很多优良的性质。比如线性，对称性(可以用在计算信号的傅里叶变换里面);
时移性：函数在时域中的时移，对应于其在频率域中附加产生的相移，而幅度频谱则保持不变;
频移性：函数在时域中乘以e^jwt，可以使整个频谱搬移w。这个也叫调制定理，通讯里面信号的频分复用需要用到这个特性(将不同的信号调制到不同的频段上同时传输);
卷积定理：时域卷积等于频域乘积;时域乘积等于频域卷积(附加一个系数)。(图像处理里面这个是个重点)
信号在频率域的表现
在频域中，频率越大说明原始信号变化速度越快;频率越小说明原始信号越平缓。当频率为0时，表示直流信号，没有变化。因此，频率的大小反应了信号的变化快慢。高频分量解释信号的突变部分，而低频分量决定信号的整体形象。
在图像处理中，频域反应了图像在空域灰度变化剧烈程度，也就是图像灰度的变化速度，也就是图像的梯度大小。对图像而言，图像的边缘部分是突变部分，变化较快，因此反应在频域上是高频分量;图像的噪声大部分情况下是高频部分;图像平缓变化部分则为低频分量。也就是说，傅立叶变换提供另外一个角度来观察图像，可以将图像从灰度分布转化到频率分布上来观察图像的特征。书面一点说就是，傅里叶变换提供了一条从空域到频率自由转换的途径。对图像处理而言，以下概念非常的重要：
图像高频分量：图像突变部分;在某些情况下指图像边缘信息，某些情况下指噪声，更多是两者的混合;
低频分量：图像变化平缓的部分，也就是图像轮廓信息
高通滤波器：让图像使低频分量抑制，高频分量通过
低通滤波器：与高通相反，让图像使高频分量抑制，低频分量通过
带通滤波器：使图像在某一部分的频率信息通过，其他过低或过高都抑制
还有个带阻滤波器，是带通的反。
模板运算与卷积定理
在时域内做模板运算，实际上就是对图像进行卷积。模板运算是图像处理一个很重要的处理过程，很多图像处理过程，比如增强/去噪(这两个分不清楚)，边缘检测中普遍用到。根据卷积定理，时域卷积等价与频域乘积。因此，在时域内对图像做模板运算就等效于在频域内对图像做滤波处理。
比如说一个均值模板，其频域响应为一个低通滤波器;在时域内对图像作均值滤波就等效于在频域内对图像用均值模板的频域响应对图像的频域响应作一个低通滤波。
图像去噪
图像去噪就是压制图像的噪音部分。因此，如果噪音是高频额，从频域的角度来看，就是需要用一个低通滤波器对图像进行处理。通过低通滤波器可以抑制图像的高频分量。但是这种情况下常常会造成边缘信息的抑制。常见的去噪模板有均值模板，高斯模板等。这两种滤波器都是在局部区域抑制图像的高频分量，模糊图像边缘的同时也抑制了噪声。还有一种非线性滤波-中值滤波器。中值滤波器对脉冲型噪声有很好的去掉。因为脉冲点都是突变的点，排序以后输出中值，那么那些最大点和最小点就可以去掉了。中值滤波对高斯噪音效果较差。
椒盐噪声：对于椒盐采用中值滤波可以很好的去除。用均值也可以取得一定的效果，但是会引起边缘的模糊。
高斯白噪声：白噪音在整个频域的都有分布，好像比较困难。
冈萨雷斯版图像处理P185：算术均值滤波器和几何均值滤波器(尤其是后者)更适合于处理高斯或者均匀的随机噪声。谐波均值滤波器更适合于处理脉冲噪声。
图像增强
有时候感觉图像增强与图像去噪是一对矛盾的过程，图像增强经常是需要增强图像的边缘，以获得更好的显示效果，这就需要增加图像的高频分量。而图像去噪是为了消除图像的噪音，也就是需要抑制高频分量。有时候这两个又是指类似的事情。比如说，消除噪音的同时图像的显示效果显著的提升了，那么，这时候就是同样的意思了。
常见的图像增强方法有对比度拉伸，直方图均衡化，图像锐化等。前面两个是在空域进行基于像素点的变换，后面一个是在频域处理。我理解的锐化就是直接在图像上加上图像高通滤波后的分量，也就是图像的边缘效果。对比度拉伸和直方图均衡化都是为了提高图像的对比度，也就是使图像看起来差异更明显一些，我想，经过这样的处理以后，图像也应该增强了图像的高频分量，使得图像的细节上差异更大。同时也引入了一些噪音。
图像傅立叶变换的物理意义
图像的频率是表征图像中灰度变化剧烈程度的指标，是灰度在平面空间上的梯度。如：大面积的沙漠在图像中是一片灰度变化缓慢的区域，对应的频率值很低;而对于地表属性变换剧烈的边缘区域在图像中是一片灰度变化剧烈的区域，对应的频率值较高。
傅立叶变换在实际中有非常明显的物理意义，设f是一个能量有限的模拟信号，则其傅立叶变换就表示f的谱。从纯粹的数学意义上看，傅立叶变换是将一个函数转换为一系列周期函数来处理的。从物理效果看，傅立叶变换是将图像从空间域转换到频率域，其逆变换是将图像从频率域转换到空间域。换句话说，傅立叶变换的物理意义是将图像的灰度分布函数变换为图像的频率分布函数，傅立叶逆变换是将图像的频率分布函数变换为灰度分布函数。
傅立叶变换以前，图像(未压缩的位图)是由对在连续空间(现实空间)上的采样得到一系列点的集合，我们习惯用一个二维矩阵表示空间上各点，则图像可由z=f(x,y)来表示。由于空间是三维的，图像是二维的，因此空间中物体在另一个维度上的关系就由梯度来表示，这样我们可以通过观察图像得知物体在三维空间中的对应关系。
为什么要提梯度?因为实际上对图像进行二维傅立叶变换得到频谱图，就是图像梯度的分布图，当然频谱图上的各点与图像上各点并不存在一一对应的关系，即使在不移频的情况下也是没有。傅立叶频谱图上我们看到的明暗不一的亮点，实际上图像上某一点与邻域点差异的强弱，即梯度的大小，也即该点的频率的大小(可以这么理解，图像中的低频部分指低梯度的点，高频部分相反)。一般来讲，梯度大则该点的亮度强，否则该点亮度弱。这样通过观察傅立叶变换后的频谱图，也叫功率图，我们首先就可以看出，图像的能量分布，如果频谱图中暗的点数更多，那么实际图像是比较柔和的(因为各点与邻域差异都不大，梯度相对较小)，反之，如果频谱图中亮的点数多，那么实际图像一定是尖锐的，边界分明且边界两边像素差异较大的。对频谱移频到原点以后，可以看出图像的频率分布是以原点为圆心，对称分布的。将频谱移频到圆心除了可以清晰地看出图像频率分布以外，还有一个好处，它可以分离出有周期性规律的干扰信号，比如正弦干扰，一副带有正弦干扰，移频到原点的频谱图上可以看出除了中心以外还存在以某一点为中心，对称分布的亮点集合，这个集合就是干扰噪音产生的，这时可以很直观的通过在该位置放置带阻滤波器消除干扰。
另外我还想说明以下几点：
1、图像经过二维傅立叶变换后，其变换系数矩阵表明：
若变换矩阵Fn原点设在中心，其频谱能量集中分布在变换系数短阵的中心附近(图中阴影区)。若所用的二维傅立叶变换矩阵Fn的原点设在左上角，那么图像信号能量将集中在系数矩阵的四个角上。这是由二维傅立叶变换本身性质决定的。同时也表明一股图像能量集中低频区域。
2 、变换之后的图像在原点平移之前四角是低频，最亮，平移之后中间部分是低频，最亮，亮度大说明低频的能量大(幅角比较大)
大致思路可能是：解析函数属调和函数一部分，而圆盘内调和函数可表示为其圆周上的限制的一个积分，这个地方好象和富变换有些联系。你这个问题很有起发，我也一直想这类问题。
拉普拉斯变换的推导途径：
1、 从数学角度：通过积分变换进行函数到函数的变换，将微分方程变为代数方程。
2、 从物理意义推导：本质上依然是将信号分解为多个正交的子信号的和(积分)，或可以从FT推广出。
从傅里叶变换导出拉普拉斯变换，可以更加清晰地解释其物理含义，并且可以将两种变换紧密地联系起来。
拉普拉斯变换提供了一种变换定义域的方法，把定义在时域上的信号(函数)映射到复频域上(要理解这句话，需要了解一下函数空间的概念--我们知道，函数定义了一种“从一个集合的元素到另一个集合的元素”的关系，而两个或以上的函数组合成的集合，就是函数空间，即函数空间也是一个集合;拉普拉斯变换的“定义域”，就是函数空间，可以说，拉普拉斯变换就是一种处理函数的函数。由于拉普拉斯变换定义得相当巧妙，所以它就具有一些奇特的特质)，而且，这是一种一一对应的关系(只要给定复频域的收敛域)，故只要给定一个时域函数(信号)，它就能通过拉普拉斯变换变换到一个复频域信号(不管这个信号是实信号还是复信号)，因而，只要我们对这个复频域信号进行处理，也就相当于对时域信号进行处理(例如设f(t)←→F(s),Re>a，则若我们对F(s)进行时延处理，得到信号F(s-z)，Re>a+Re,那么就相当于我们给时域函数乘以一个旋转因子e^zt，即f(t)e^zt←→F(s-z),Re>a+Re;只要对F(s-z)进行反变换，就可以得到f(t)e^zt)。 拉普拉斯变换被用于求解微分方程，主要是应用拉普拉斯变换的几个性质，使求解微分方程转变为求解代数方程(因为求解代数方程总比求解微分方程容易得多!而且，(可以很方便地)对求解结果进行拉普拉斯反变换从而得到原微分方程的解)。
我们总可以容易地画出实变函数的图像(绝大多数函数的确如此)，但我们难以画出一个复变函数的图象，这也许是拉普拉斯变换比较抽象的原因之一;而另外一个原因，就是拉普拉斯变换中的复频率s没有明确的物理意义。
关于特征根和复数，建议提问者再去看看书中的定义，应该不难理解。
-----------------------------------------------------------------------------------------------------------
示例图像用32*32像素，16色灰度
问题补充：
void Fft(char d1[],char f[],char fi[]) { int i,j,x,y; char t; for(i=0;i<32;i++) { for(j=0;j<32;i++) { t=pow(-1,i+j)*d1; for(x=0;x<32;i++) { for(y=0;y<32;i++) { f+=t*cos((2*pi*i*x+2*pi*j*y)/32); fi+=t*sin((2*pi*i*x+2*pi*j*y)/32); } } d1=sqrt(pow(f,2)+pow(fi,2)); putpixel(j,i,d1); } } } void lowpass(char f[],char fi[]) { int i,j,dx; int d0=5; for(i=0;i<32;i++) { for(j=0;j<32;i++) { dx=sqrt(pow(i-16,2)+pow(j-16,2)); if(dx>d0) { f=0; fi=0; } } } }
转自：http://blog.sina.com.cn/s/blog_79496d6b0100rra1.html
快速傅立叶变换的意义及应用
为什么要进行傅里叶变换，其物理意义是什么?傅立叶变换是数字信号处理领域一种很重要的算法。要知道傅立叶变换算法的意义，首先要了解傅立叶原理的意义。傅立叶原理表明：任何连续测量的时序或信号，都可以表示为不同频率的正弦波信号的无限叠加。而根据该原理创立的傅立叶变换算法利用直接测量到的原始信号，以累加方式来计算该信号中不同正弦波信号的频率、振幅和相位。
和傅立叶变换算法对应的是反傅立叶变换算法。该反变换从本质上说也是一种累加处理，这样就可以将单独改变的正弦波信号转换成一个信号。
因此，可以说，傅立叶变换将原来难以处理的时域信号转换成了易于分析的频域信号(信号的频谱)，可以利用一些工具对这些频域信号进行处理、加工。最后还可以利用傅立叶反变换将这些频域信号转换成时域信号。
从现代数学的眼光来看，傅里叶变换是一种特殊的积分变换。它能将满足一定条件的某个函数表示成正弦基函数的线性组合或者积分。在不同的研究领域，傅里叶变换具有多种不同的变体形式，如连续傅里叶变换和离散傅里叶变换。
在数学领域，尽管最初傅立叶分析是作为热过程的解析分析的工具，但是其思想方法仍然具有典型的还原论和分析主义的特征。"任意"的函数通过一定的分解，都能够表示为正弦函数的线性组合的形式，而正弦函数在物理上是被充分研究而相对简单的函数类：1. 傅立叶变换是线性算子,若赋予适当的范数,它还是酉算子;2. 傅立叶变换的逆变换容易求出,而且形式与正变换非常类似;3. 正弦基函数是微分运算的本征函数,从而使得线性微分方程的求解可以转化为常系数的代数方程的求解.在线性时不变杂的卷积运算为简单的乘积运算,从而提供了计算卷积的一种简单手段;5. 离散形式的傅立叶的物理系统内,频率是个不变的性质,从而系统对于复杂激励的响应可以通过组合其对不同频率正弦信号的响应来获取;4. 著名的卷积定理指出:傅立叶变换可以化复变换可以利用数字计算机快速的算出(其算法称为快速傅立叶变换算法(FFT))。
正是由于上述的良好性质,傅里叶变换在物理学、数论、组合数学、信号处理、概率、统计、密码学、声学、光学等领域都有着广泛的应用。
2、图像傅立叶变换的物理意义
图像的频率是表征图像中灰度变化剧烈程度的指标，是灰度在平面空间上的梯度。如：大面积的沙漠在图像中是一片灰度变化缓慢的区域，对应的频率值很低;而对于地表属性变换剧烈的边缘区域在图像中是一片灰度变化剧烈的区域，对应的频率值较高。傅立叶变换在实际中有非常明显的物理意义，设f是一个能量有限的模拟信号，则其傅立叶变换就表示f的谱。从纯粹的数学意义上看，傅立叶变换是将一个函数转换为一系列周期函数来处理的。从物理效果看，傅立叶变换是将图像从空间域转换到频率域，其逆变换是将图像从频率域转换到空间域。换句话说，傅立叶变换的物理意义是将图像的灰度分布函数变换为图像的频率分布函数，傅立叶逆变换是将图像的频率分布函数变换为灰度分布函数
傅立叶变换以前，图像(未压缩的位图)是由对在连续空间(现实空间)上的采样得到一系列点的集合，我们习惯用一个二维矩阵表示空间上各点，则图像可由z=f(x,y)来表示。由于空间是三维的，图像是二维的，因此空间中物体在另一个维度上的关系就由梯度来表示，这样我们可以通过观察图像得知物体在三维空间中的对应关系。为什么要提梯度?因为实际上对图像进行二维傅立叶变换得到频谱图，就是图像梯度的分布图，当然频谱图上的各点与图像上各点并不存在一一对应的关系，即使在不移频的情况下也是没有。傅立叶频谱图上我们看到的明暗不一的亮点，实际上图像上某一点与邻域点差异的强弱，即梯度的大小，也即该点的频率的大小(可以这么理解，图像中的低频部分指低梯度的点，高频部分相反)。一般来讲，梯度大则该点的亮度强，否则该点亮度弱。这样通过观察傅立叶变换后的频谱图，也叫功率图，我们首先就可以看出，图像的能量分布，如果频谱图中暗的点数更多，那么实际图像是比较柔和的(因为各点与邻域差异都不大，梯度相对较小)，反之，如果频谱图中亮的点数多，那么实际图像一定是尖锐的，边界分明且边界两边像素差异较大的。对频谱移频到原点以后，可以看出图像的频率分布是以原点为圆心，对称分布的。将频谱移频到圆心除了可以清晰地看出图像频率分布以外，还有一个好处，它可以分离出有周期性规律的干扰信号，比如正弦干扰，一副带有正弦干扰，移频到原点的频谱图上可以看出除了中心以外还存在以某一点为中心，对称分布的亮点集合，这个集合就是干扰噪音产生的，这时可以很直观的通过在该位置放置带阻滤波器消除干扰
快速傅氏变换 英文名是fast Fourier transform
快速傅氏变换(FFT)是离散傅氏变换(DFT)的快速算法，它是根据离散傅氏变换的奇、偶、虚、实等特性，对离散傅立叶变换的算法进行改进获得的。它对傅氏变换的理论并没有新的发现，但是对于在计算机系统或者说数字系统中应用离散傅立叶变换，可以说是进了一大步。
设x(n)为N项的复数序列，由DFT变换，任一X(m)的计算都需要N次复数乘法和N-1次复数加法，而一次复数乘法等于四次实数乘法和两次实数加法，一次复数加法等于两次实数加法，即使把一次复数乘法和一次复数加法定义成一次“运算”(四次实数乘法和四次实数加法)，那么求出N项复数序列的X(m),即N点DFT变换大约就需要N2次运算。当N=1024点甚至更多的时候，需要N2=1048576次运算，在FFT中，利用WN的周期性和对称性，把一个N项序列(设N=2k,k为正整数)，分为两个N/2项的子序列，每个N/2点DFT变换需要(N/2)2次运算，再用N次运算把两个N/2点的DFT变换组合成一个N点的DFT变换。这样变换以后，总的运算次数就变成N+2(N/2)2=N+N2/2。继续上面的例子，N=1024时，总的运算次数就变成了525312次，节省了大约50%的运算量。而如果我们将这种“一分为二”的思想不断进行下去，直到分成两两一组的DFT运算单元，那么N点的DFT变换就只需要Nlog2N次的运算，N在1024点时，运算量仅有10240次，是先前的直接算法的1%，点数越多，运算量的节约就越大，这就是FFT的优越性。
小波分析 (Wavelet)
小波分析是当前数学中一个迅速发展的新领域，它同时具有理论深刻和应用十分广泛的双重意义。
小波变换的概念是由法国从事石油信号处理的工程师J.Morlet在1974年首先提出的，通过物理的直观和信号处理的实际需要经验的建立了反演公式，当时未能得到数学家的认可。正如1807年法国的热学工程师J.B.J.Fourier提出任一函数都能展开成三角函数的无穷级数的创新概念未能得到著名数学家J.L.Lagrange，P.S.Laplace以及A.M.Legendre的认可一样。幸运的是，早在七十年代，A.Calderon表示定理的发现、Hardy空间的原子分解和无条件基的深入研究为小波变换的诞生做了理论上的准备，而且J.O.Stromberg还构造了历史上非常类似于现在的小波基;1986年著名数学家Y.Meyer偶然构造出一个真正的小波基，并与S.Mallat合作建立了构造小波基的同意方法枣多尺度分析之后，小波分析才开始蓬勃发展起来，其中比利时女数学家I.Daubechies撰写的《小波十讲(Ten Lectures on Wavelets)》对小波的普及起了重要的推动作用。它与Fourier变换、窗口Fourier变换(Gabor变换)相比，这是一个时间和频率的局域变换，因而能有效的从信号中提取信息，通过伸缩和平移等运算功能对函数或信号进行多尺度细化分析(Multiscale Analysis)，解决了Fourier变换不能解决的许多困难问题，从而小波变化被誉为“数学显微镜”，它是调和分析发展史上里程碑式的进展。
小波(Wavelet)这一术语，顾名思义，“小波”就是小的波形。所谓“小”是指它具有衰减性;而称之为“波”则是指它的波动性，其振幅正负相间的震荡形式。与Fourier变换相比，小波变换是时间(空间)频率的局部化分析，它通过伸缩平移运算对信号(函数)逐步进行多尺度细化，最终达到高频处时间细分，低频处频率细分，能自动适应时频信号分析的要求，从而可聚焦到信号的任意细节，解决了Fourier变换的困难问题，成为继Fourier变换以来在科学方法上的重大突破。有人把小波变换称为“数学显微镜”。
小波分析的应用是与小波分析的理论研究紧密地结合在一起地。现在，它已经在科技信息产业领域取得了令人瞩目的成就。电子信息技术是六大高新技术中重要的一个领域，它的重要方面是图像和信号处理。现今，信号处理已经成为当代科学技术工作的重要部分，信号处理的目的就是：准确的分析、诊断、编码压缩和量化、快速传递或存储、精确地重构(或恢复)。从数学地角度来看，信号与图像处理可以统一看作是信号处理(图像可以看作是二维信号)，在小波分析地许多分析的许多应用中，都可以归结为信号处理问题。现在，对于其性质随实践是稳定不变的信号，处理的理想工具仍然是傅立叶分析。但是在实际应用中的绝大多数信号是非稳定的，而特别适用于非稳定信号的工具就是小波分析。
小波分析是当前应用数学和工程学科中一个迅速发展的新领域，经过近10年的探索研究，重要的数学形式化体系已经建立，理论基础更加扎实。与Fourier变换相比，小波变换是空间(时间)和频率的局部变换，因而能有效地从信号中提取信息。通过伸缩和平移等运算功能可对函数或信号进行多尺度的细化分析，解决了Fourier变换不能解决的许多困难问题。小波变换联系了应用数学、物理学、计算机科学、信号与信息处理、图像处理、地震勘探等多个学科。数学家认为，小波分析是一个新的数学分支，它是泛函分析、Fourier分析、样调分析、数值分析的完美结晶;信号和信息处理专家认为，小波分析是时间—尺度分析和多分辨分析的一种新技术，它在信号分析、语音合成、图像识别、计算机视觉、数据压缩、地震勘探、大气与海洋波分析等方面的研究都取得了有科学意义和应用价值的成果。
事实上小波分析的应用领域十分广泛，它包括：数学领域的许多学科;信号分析、图像处理;量子力学、理论物理;军事电子对抗与武器的智能化;计算机分类与识别;音乐与语言的人工合成;医学成像与诊断;地震勘探数据处理;大型机械的故障诊断等方面;例如，在数学方面，它已用于数值分析、构造快速数值方法、曲线曲面构造、微分方程求解、控制论等。在信号分析方面的滤波、去噪声、压缩、传递等。在图像处理方面的图像压缩、分类、识别与诊断，去污等。在医学成像方面的减少B超、CT、核磁共振成像的时间，提高分辨率等。
(1)小波分析用于信号与图像压缩是小波分析应用的一个重要方面。它的特点是压缩比高，压缩速度快，压缩后能保持信号与图像的特征不变，且在传递中可以抗干扰。基于小波分析的压缩方法很多，比较成功的有小波包最好基方法，小波域纹理模型方法，小波变换零树压缩，小波变换向量压缩等。
(2)小波在信号分析中的应用也十分广泛。它可以用于边界的处理与滤波、时频分析、信噪分离与提取弱信号、求分形指数、信号的识别与诊断以及多尺度边缘检测等。
(3)在工程技术等方面的应用。包括计算机视觉、计算机图形学、曲线设计、湍流、远程宇宙的研究与生物医学方面。
傅立叶变换的基函数是sin cos函数，也就是说它只是用这2个去逼近原信号，对于原信号不一定都与正弦和余弦函数相似，因此小波基函数的多样性能使逼近的权重因子更加准确.正是因为FFT的频域局部性不强，而且时频分开，所以才有了以后的方法――当然包括了小波!
转自：http://blog.sina.com.cn/s/blog_62dc34180100fs33.html
功率谱密度相关方法的MATLAB实现
1. 基本方法周期图法是直接将信号的采样数据x(n)进行Fourier变换求取功率谱密度估计的方法。假定有限长随机信号序列为x(n)。它的Fourier变换和功率谱密度估计存在下面的关系：
式中，N为随机信号序列x(n)的长度。在离散的频率点f=kΔf,有：
其中，FFT为对序列x(n)的Fourier变换，由于FFT的周期为N，求得的功率谱估计以N为周期，因此这种方法称为周期图法。下面用例子说明如何采用这种方法进行功率谱
用有限长样本序列的Fourier变换来表示随机序列的功率谱，只是一种估计或近似，不可避免存在误差。为了减少误差，使功率谱估计更加平滑，可采用分段平均周期图法(Bartlett法)、加窗平均周期图法(Welch法)等方法加以改进。
2. 分段平均周期图法(Bartlett法)
将信号序列x(n)，n=0,1,…,N-1，分成互不重叠的P个小段，每小段由m个采样值，则P*m=N。对每个小段信号序列进行功率谱估计，然后再取平均作为整个序列x(n)的功率谱估计。
平均周期图法还可以对信号x(n)进行重叠分段，如按2:1重叠分段，即前一段信号和后一段信号有一半是重叠的。对每一小段信号序列进行功率谱估计，然后再取平均值作为整个序列x(n)的功率谱估计。这两种方法都称为平均周期图法，一般后者比前者好。程序运行结果为图9-5，上图采用不重叠分段法的功率谱估计，下图为2:1重叠分段的功率谱估计，可见后者估计曲线较为平滑。与上例比较，平均周期图法功率谱估计具有明显效果(涨落曲线靠近0dB)。
3.加窗平均周期图法
加窗平均周期图法是对分段平均周期图法的改进。在信号序列x(n)分段后，用非矩形窗口对每一小段信号序列进行预处理，再采用前述分段平均周期图法进行整个信号序列x(n)的功率谱估计。由窗函数的基本知识(第7章)可知，采用合适的非矩形窗口对信号进行处理可减小“频谱泄露”，同时可增加频峰的宽度，从而提高频谱分辨率。
其中上图采用无重叠数据分段的加窗平均周期图法进行功率谱估计，而下图采用重叠数据分段的加窗平均周期图法进行功率谱估计，显然后者是更佳的，信号谱峰加宽，而噪声谱均在0dB附近，更为平坦(注意采用无重叠数据分段噪声的最大的下降分贝数大于5dB，而重叠数据分段周期图法噪声的最大下降分贝数小于5dB)。
4. Welch法估计及其MATLAB函数
Welch功率谱密度就是用改进的平均周期图法来求取随机信号的功率谱密度估计的。Welch法采用信号重叠分段、加窗函数和FFT算法等计算一个信号序列的自功率谱估计(PSD如上例中的下半部分的求法)和两个信号序列的互功率谱估计(CSD)。
MATLAB信号处理工具箱函数提供了专门的函数PSD和CSD自动实现Welch法估计，而不需要自己编程。
(1) 函数psd利用Welch法估计一个信号自功率谱密度，函数调用格式为：
]=psd(x[,Nfft,Fs,window,Noverlap,’dflag’])
式中，x为信号序列;Nfft为采用的FFT长度。这一值决定了功率谱估计速度，当Nfft采用2的幂时，程序采用快速算法;Fs为采样频率;Window定义窗函数和x分段序列的长度。窗函数长度必须小于或等于Nfft，否则会给出错误信息;Noverlap为分段序列重叠的采样点数(长度)，它应小于Nfft;dflag为去除信号趋势分量的选择项：’linear’,去除线性趋势分量，’mean’去除均值分量，’none’不做去除趋势处理。Pxx为信号x的自功率谱密度估计。f为返回的频率向量，它和Pxx对应，并且有相同长度。
在psd函数调用格式中，缺省值为：Nfft=min(256,length(x)),Fs=2Hz, window=hanning(Nfft),noverlap=0. 若x是实序列，函数psd仅计算频率为正的功率
注意程序前半部分中频率向量f的创建方法。它与函数psd的输出Pxx长度的关系如下：若x为实序列，当Nfft为奇数时，f=(0:(Nfft+1)/2-1)/Nfft;当Nfft为偶数时，f=(0:Nfft/2)/Nfft。
函数还有一种缺省返回值的调用格式，用于直接绘制信号序列x的功率谱估计曲线。
函数还可以计算带有置信区间的功率谱估计，调用格式为：
=psd(x,Nfft,Fs,window,Noverlap,p)
式中，p为置信区间，0<=p<=1。
由此可知，滤波器输入白噪声序列的输出信号的功率谱或自相关可以确定滤波器的频率特性。
(2)函数csd利用welch法估计两个信号的互功率谱密度，函数调用格式为：
]=csd(x,y,Nfft,Fs,window,Noverlap,’dflag’)
]=csd(x,y,Nfft,Fs,window,Noverlap,p)
这里，x,y为两个信号序列;Pxy为x,y的互功率谱估计;其他参数的意义同自功率谱函数psd。
可以看到，两个白噪声信号的互功率谱(上图)杂乱无章，看不出周期成分，大部分功率谱在-5dB以下。然而白噪声与带有噪声的周期信号的功率谱在其周期(频率为1000Hz)处有一峰值，清楚地表明了周期信号的周期或频率。因此，利用未知信号与白噪声信号的互功率谱也可以检测未知信号中所含有的频率成分。
5 多 窗 口 法
多窗口法(Multitaper method,简称MTM法)利用多个正交窗口(Tapers)获得各自独立的近似功率谱估计，然后综合这些估计得到一个序列的功率谱估计。相对于普通的周期图法，这种功率谱估计具有更大的自由度，并在估计精度和估计波动方面均有较好的效果。普通的功率谱估计只利用单一窗口，因此在序列始端和末端均会丢失相关信息，而且无法找回。而MTM法估计增加窗口使得丢失的信息尽量减少。
MTM法简单地采用一个参数：时间带宽积(Time-bandwidth product)NW，这个参数用以定义计算功率谱所用窗的数目，为2*NW-1。NW越大，功率谱计算次数越多，时间域分辨率越高，而频率域分辨率降低，使得功率谱估计的波动减小。随着NW增大，每次估计中谱泄漏增多，总功率谱估计的偏差增大。对于每一个数据组，通常有一个最优的NW使得在估计偏差和估计波动两方面求得折中，这需要在程序中反复调试来获得。
MATLAB信号处理工具箱中函数PMTM就是采用MTM法估计功率谱密度。函数调用格式为：
]=pmtm(x[,nw,Nfft,Fs])
式中，x为信号序列;nw为时间带宽积，缺省值为4。通常可取2，5/2，3，7/2;Nfft为FFT长度;Fs为采样频率。
上面的函数还可以通过无返回值而绘出置信区间，如pmtm(x,nw,Nfft,Fs,’option’,p)绘制带置信区间的功率谱密度估计曲线，0<=p<=1。
6 最 大 熵 法(Maxmum entropy method， MEM法)
如上所述，周期图法功率谱估计需要对信号序列“截断”或加窗处理，其结果是使估计的功率谱密度为信号序列真实谱和窗谱的卷积，导致误差的产生。
最大熵功率谱估计的目的是最大限度地保留截断后丢失的“窗口”以外信号的信息，使估计谱的熵最大。主要方法是以已知的自相关序列rxx(0),rxx(1),…,rxx(p)为基础，外推自相关序列rxx(p+1),rxx(p+2),…，保证信息熵最大。
最大熵功率谱估计法假定随机过程是平稳高斯过程，可以证明，随机信号的最大熵谱与AR自回归(全极点滤波器)模型谱是等价的。
MATLAB信号处理工具箱提供最大熵功率谱估计函数pmem，其调用格式为：
=pmem(x,p,Nfft,Fs,’xcorr’)
式中，x为输入信号序列或输入相关矩阵;p为全极点滤波器阶次;a为全极点滤波器模型系数向量;’xcorr’是把x认为是相关矩阵。
比较最大熵功率谱估计(MEM)和改进的平均周期图功率谱估计，可见，MEM法估计的功率谱曲线较光滑。在这一方法中，MEM法选定全极点滤波器的阶数取得越大，能够获得的窗口外的信息越多，但计算量也越大，需要根据情况折中考虑。
7 多信号分类法
MATLAB信号处理工具箱还提供另一种功率谱估计函数pmusic。该函数执行多信号分类法(multiple signal classification， Music法)。将数据自相关矩阵看成由信号自相关矩阵和噪声自相关矩阵两部分组成，即数据自相关矩阵R包含有两个子空间信息：信号子空间和噪声子空间。这样，矩阵特征值向量(Eigen vector)也可分为两个子空间：信号子空间和噪声子空间。为了求得功率谱估计，函数pmusic计算信号子空间和噪声子空间的特征值向量函数，使得在周期信号频率处函数值最大，功率谱估计出现峰值，而在其他频率处函数值最小。其调用格式为：
=pmusic(x,p[[,thresh],Nfft,Fs,window,Noverlap])
式中,x为输入信号的向量或矩阵;p为信号子空间维数;thresh为阈值，其他参数的意义与函数psd相同。
功率谱密度相关方法的MATLAB实现
%%
%分段平均周期图法(Bartlett法)
%运用信号不重叠分段估计功率谱
Nsec=256;n=0:sigLength-1;t=n/Fs; %数据点数，分段间隔，时间序列
pxx1=abs(fft(y(1:256),Nsec).^2)/Nsec; %第一段功率谱
pxx2=abs(fft(y(257:512),Nsec).^2)/Nsec; %第二段功率谱
pxx3=abs(fft(y(515:768),Nsec).^2)/Nsec; %第三段功率谱
pxx4=abs(fft(y(769:1024),Nsec).^2)/Nsec; %第四段功率谱
Pxx=10*log10((pxx1+pxx2+pxx3+pxx4)/4); %平均得到整个序列功率谱
f=(0:length(Pxx)-1)*Fs/length(Pxx); %给出功率谱对应的频率
%%plot(f(1:Nsec/2),Pxx(1:Nsec/2)); %绘制功率谱曲线
figure,plot(f(1:Nsec),Pxx(1:Nsec)); %绘制功率谱曲线
xlabel('频率/Hz');ylabel('功率谱 /dB');
title('平均周期图(无重叠) N=4*256');
grid on
%%
%运用信号重叠分段估计功率谱
pxx1=abs(fft(y(1:256),Nsec).^2)/Nsec; %第一段功率谱
pxx2=abs(fft(y(129:384),Nsec).^2)/Nsec; %第二段功率谱
pxx3=abs(fft(y(257:512),Nsec).^2)/Nsec; %第三段功率谱
pxx4=abs(fft(y(385:640),Nsec).^2)/Nsec; %第四段功率谱
pxx5=abs(fft(y(513:768),Nsec).^2)/Nsec; %第五段功率谱
pxx6=abs(fft(y(641:896),Nsec).^2)/Nsec; %第六段功率谱
pxx7=abs(fft(y(769:1024),Nsec).^2)/Nsec; %第七段功率谱
Pxx=10*log10((pxx1+pxx2+pxx3+pxx4+pxx5+pxx6+pxx7)/7); %功率谱平均并转化为dB
f=(0:length(Pxx)-1)*Fs/length(Pxx); %频率序列
figure,plot(f(1:Nsec/2),Pxx(1:Nsec/2)); %绘制功率谱曲线
xlabel('频率/Hz'); ylabel('功率谱/dB');
title('平均周期图(重叠一半) N=1024');
grid on
%%
Nsec=256;n=0:sigLength-1;t=n/Fs; %数据长度，分段数据长度、时间序列
w=hanning(256); %采用的窗口数据
%采用不重叠加窗方法的功率谱估计
pxx1=abs(fft(w.*y(1:256),Nsec).^2)/norm(w)^2; %第一段加窗振幅谱平方
pxx2=abs(fft(w.*y(257:512),Nsec).^2)/norm(w)^2; %第二段加窗振幅谱平方
pxx3=abs(fft(w.*y(513:768),Nsec).^2)/norm(w)^2; %第三段加窗振幅谱平方
pxx4=abs(fft(w.*y(769:1024),Nsec).^2)/norm(w)^2; %第四段加窗振幅谱平方
Pxx=10*log10((pxx1+pxx2+pxx3+pxx4)/4); %求得平均功率谱，转换为dB
f=(0:length(Pxx)-1)*Fs/length(Pxx); %求得频率序列
figure
subplot(2,1,1),plot(f(1:Nsec/2),Pxx(1:Nsec/2)); %绘制功率谱曲线
xlabel('频率/Hz');ylabel('功率谱/dB');
title('加窗平均周期图(无重叠) N=4*256');
grid on
%采用重叠加窗方法的功率谱估计
pxx1=abs(fft(w.*y(1:256),Nsec).^2)/norm(w)^2; %第一段加窗振幅谱平方
pxx2=abs(fft(w.*y(129:384),Nsec).^2)/norm(w)^2; %第二段加窗振幅谱平方
pxx3=abs(fft(w.*y(257:512),Nsec).^2)/norm(w)^2; %第三段加窗振幅谱平方
pxx4=abs(fft(w.*y(385:640),Nsec).^2)/norm(w)^2; %第四段加窗振幅谱平方
pxx5=abs(fft(w.*y(513:768),Nsec).^2)/norm(w)^2; %第五段加窗振幅谱平方
pxx6=abs(fft(w.*y(641:896),Nsec).^2)/norm(w)^2; %第六段加窗振幅谱平方
pxx7=abs(fft(w.*y(769:1024),Nsec).^2)/norm(w)^2; %第七段加窗振幅谱平方
Pxx=10*log10((pxx1+pxx2+pxx3+pxx4+pxx5+pxx6+pxx7)/7);%平均功率谱转换为dB
f=(0:length(Pxx)-1)*Fs/length(Pxx); %频率序列
subplot(2,1,2),plot(f(1:Nsec/2),Pxx(1:Nsec/2)); %绘制功率谱曲线
xlabel('频率/Hz');ylabel('功率谱/dB');
title('加窗平均周期图(重叠一半)N=1024');
grid on
%%
%4分段平均周期图法(hanning窗)
Nsec=256;
n=0:sigLength-1;
t=n/Fs;
w=hanning(256);
Pxx1=abs(fft(w.*y(1:256),Nsec).^2)/Nsec;
Pxx2=abs(fft(w.*y(257:512),Nsec).^2)/Nsec;
Pxx3=abs(fft(w.*y(513:768),Nsec).^2)/Nsec;
Pxx4=abs(fft(w.*y(769:1024),Nsec).^2)/Nsec;
Pxx=10*log10((Pxx1+Pxx2+Pxx3+Pxx4)/4);
f=(0:length(Pxx)-1)*Fs/length(Pxx);
figure
subplot(2,1,1)
plot(f,Pxx);
xlabel('频率/Hz');ylabel('功率谱/dB');
title('Averaged Modified Periodogram (none overlap) N=4*256');
grid
%4分段(2：1重叠)平均周期图法(hanning窗)
Nsec=256;
n=0:sigLength-1;
t=n/Fs;
w=hanning(256);
Pxx1=abs(fft(w.*y(1:256),Nsec).^2)/Nsec;
Pxx2=abs(fft(w.*y(129:384),Nsec).^2)/Nsec;
Pxx3=abs(fft(w.*y(257:512),Nsec).^2)/Nsec;
Pxx4=abs(fft(w.*y(385:640),Nsec).^2)/Nsec;
Pxx5=abs(fft(w.*y(513:768),Nsec).^2)/Nsec;
Pxx6=abs(fft(w.*y(641:896),Nsec).^2)/Nsec;
Pxx7=abs(fft(w.*y(769:1024),Nsec).^2)/Nsec;
Pxx=10*log10((Pxx1+Pxx2+Pxx3+Pxx4+Pxx5+Pxx6+Pxx7)/7);
f=(0:length(Pxx)-1)*Fs/length(Pxx);
subplot(2,1,2);plot(f,Pxx);
xlabel('频率/Hz');ylabel('Power Spectrum (dB)');
title('Averaged Modified Periodogram (half overlap) N=1024');
grid
%%
%PSD_WELCH方法
%采样频率
Nfft=256;n=0:sigLength-1;t=n/Fs; %数据长度、时间序列
window=hanning(256); %选用的窗口
noverlap=128; %分段序列重叠的采样点数(长度)
dflag='none'; %不做趋势处理
=psd(y,Nfft,Fs,window,noverlap,0.95); %功率谱估计,并以0.95的置信度给出置信区间，无返回值是绘制出置信区间
figure
plot(f,10*log10(Pxx)); %绘制功率谱
xlabel('频率/Hz');ylabel('功率谱/dB');
title('PSD—Welch方法'); grid on
%%
%最大熵法(MEM法)
Nfft=256;n=0:sigLength-1;t=n/Fs; %数据长度、分段长度和时间序列
window=hanning(256); %采用窗口
=pmem(x,20,Nfft,Fs); %采用最大熵法，采用滤波器阶数14，估计功率谱
figure,subplot(2,1,1),plot(f,10*log10(Pxx1)); %绘制功率谱
xlabel('频率/Hz');ylabel('功率谱/dB');
title('最大熵法 Order=20原始信号功率谱');grid on
=pmem(y0,20,Nfft,Fs); %采用最大熵法，采用滤波器阶数14，估计功率谱
subplot(2,1,2),plot(f,10*log10(Pxx1)); %绘制功率谱
xlabel('频率/Hz');ylabel('功率谱/dB');axis()
title('最大熵法 Order=20滤波后的信号功率谱');grid on
%%
%%功率谱密度
%PSD_WELCH方法
Nfft=512;n=0:sigLength-1;t=n/Fs; %数据长度、时间序列
window=hanning(256); %选用的窗口
noverlap=128; %分段序列重叠的采样点数(长度)
dflag='none'; %不做趋势处理
=psd(x,Nfft,Fs,window,noverlap,0.95); %功率谱估计,并以0.95的置信度给出置信区间，无返回值是绘制出置信区间
figure;subplot(211);
plot(f,10*log10(Pxx)); %绘制功率谱
xlabel('频率/Hz');ylabel('功率谱/dB');
grid on;title('PSD—Welch方法的原始信号功率谱')
subplot(212)
=psd(y0,Nfft,Fs,window,noverlap,0.95); %功率谱估计,并以0.95的置信度给出置信区间，无返回值是绘制出置信区间
plot(f,10*log10(Pxx)); %绘制功率谱
xlabel('频率/Hz');ylabel('功率谱/dB');axis()
grid on;title('PSD—Welch方法的滤波后的信号功率谱')
%%
%用多窗口法(MTM)
n=0:sigLength-1;t=n/Fs; %数据长度、分段数据长度，时间序列
=pmtm(x,2,Nfft,Fs); %用多窗口法(NW=4)估计功率谱
figure;subplot(2,1,1),plot(f,10*log10(Pxx1)); %绘制功率谱
xlabel('频率/Hz');ylabel('功率谱/dB');
title('多窗口法(MTM) nw=2原始信号功率谱');grid on
=pmtm(y0,2,Nfft,Fs); %用多窗口法(NW=2)估计功率谱
subplot(2,1,2),plot(f,10*log10(Pxx)); %绘制功率谱
xlabel('频率/Hz');ylabel('功率谱/dB');axis()
title('多窗口法(MTM) nw=2滤波后的信号功率谱');grid on
%%
%采用Welch方法估计功率谱
noverlap=128; %重叠数据
dflag='none'; Nfft=1024;
figure;subplot(2,1,1)
psd(x,Nfft,Fs,window,noverlap,dflag); %采用Welch方法估计功率谱
xlabel('频率/Hz');ylabel('功率谱/dB')
title('Welch方法原始信号功率谱');grid on
subplot(2,1,2)
psd(y0,Nfft,Fs,window,noverlap,dflag); %采用Welch方法估计功率谱
xlabel('频率/Hz');ylabel('功率谱/dB'),axis()
title('Welch方法滤波后的信号功率谱');grid on
pmusic(y0,,Nfft,Fs,32,16); %采用多信号分类法估计功率谱
xlabel('频率/Hz'); ylabel('功率谱/dB')
title('通过MUSIC法估计的伪谱')
转自：http://blog.sina.com.cn/s/blog_79ecf6980100vcqn.html
一个FFT例子解析,利用CImg进行快速傅立叶变换
// Drawing Demo
// All Rights Reserved by Woo
#include "stdafx.h"
#include "D:\CImg\CImg.h"
#include
using namespace std;
using namespace cimg_library;
// The lines below are necessary when using a
// non-standard compiler such as visualcpp6.
#ifdef cimg_use_visualcpp6
#define std
#endif
#ifdef min
#undef min
#undef max
#endif
// Define circumference ratio constant
#define PI 3.14159265358
int main(int argc, char* argv[])
{
// Frequency
const int feq = 25;
// Sampling rate
const int sample_rate = 1024;
// Ampitude
const float ampitude = 50.0f;
// Colors
const unsigned char red[] = { 255, 0, 0 }, yellow[] = { 255, 255, 0 },
green[] = { 0, 255, 0 };
// 1-D signal
CImg signal(1, 1024, 1, 1, 0);
// Filling signal with values
cimg_forXY( signal, x, y )
{
signal( x, y ) = ampitude*cos( 2*PI*feq*( y/ 1023.0 ) ) + 125;
}
// Create a black image as drawing panel
CImg signal_visu( 800, 300, 1, 3, 0 );
// Draw grid and axis in 'drawing panel'
signal_visu.draw_grid( 50, 50, 0, 0, 0, 0, yellow, 0.5 ).draw_axis( 0, 1000, 200, -200, green, 1.0f );
// Draw the signal with cubic spline
signal_visu.draw_graph( signal, red, 1, 2, 0, 255, 0 );
// Display our drawing result
CImgDisplay main_disp( signal_visu, "signal" );
// FFT result container
CImgList fft = signal.get_FFT( 'y' );
// Create a drawing panel
CImg fft_visu( 800, 300, 1, 3, 0 );
// Draw grid and axis in 'drawing panel'
fft_visu.draw_grid( 50, 50, 0, 0, 0, 0, yellow, 0.5 ).draw_axis( 0, 1000, 300, 0, green, 1.0f );;
// Power of FFT. I scale the power in order to display it properly
CImg power = ( fft(0).get_pow(2.0) + fft(1).get_pow(2.0) ).sqrt() * 0.05;
// Draw the signal with cubic spline
fft_visu.draw_graph( power , red, 1, 2, 0, 255, 0 );
CImgDisplay fft_disp( fft_visu, "FFT: Magitude" );
CImg phase(1, 1024, 1, 1, 0.0f);
cimg_forXY( phase, x1, y1 )
{
phase(x1, y1) = atan(abs((fft(0))(x1, y1) / (fft(1))(x1, y1)));
}
CImg phase_visu( 800, 300, 1, 3, 0 );
// Draw grid and axis in 'drawing panel'
phase_visu.draw_grid( 50, 50, 0, 0, 0, 0, yellow, 0.5 ).draw_axis( 0, 1000, 300, 0, green, 1.0f );;
// Draw the signal with cubic spline
phase_visu.draw_graph( phase , red, 1, 2, 0, 255, 0 );
CImgDisplay phase_disp( phase_visu, "FFT: Phase");
while ( !main_disp.is_closed )
{
main_disp.wait();
}
return 0;
}
今天给大家介绍的例子是利用CImg进行快速傅立叶变换。傅立叶变换被广泛的应用于信号处理行业，为信号分析提供了一种强有力的分析方法。老傅在他的一篇文章中，提出周期性的函数都可以利用正弦函数来表示。当时，审稿人中有两位当时大名顶顶的数学家：拉格朗日与拉普拉斯，这两个人对于此文章是否可以发表持相反意见。拉格朗日支持方认为该文章不能发表，他指出该方法不能表示带有边角(corner)的函数，而拉普拉斯和其他学者则支持发表。最后，迫于拉格朗日的权威(牛人就是牛人)，这篇文章在他西去之后的15年才得以发表。不过这段时间，傅立叶也没有闲着，也就没有太在意这件事情。是的，人有的时候不能太在意一些东西，随遇而安最好。到底上面两位牛人谁的观点对呢?这要一分为二看：站在老拉的角度来看，这个方法是有‘问题’，但是我们可以无限逼近元函数，使得它们之间的能量几乎为零，这样一来，傅立叶变换是可行的。此外，'对于离散信号来说，傅立叶信号是准确的，就像7 = 3 + 4一样'.(1)
简单的说，傅立叶变换其实就是统计某一频率的信号在待分析信号中出现的频率(次数)，貌似和概率分布差不多(个人观点)。傅立叶变换可以通过多种方法来解决，例如解线性方程组，求得每一幅值大小。这种方法简单直接，但是耗时巨大。我们平时所要分析的数据量往往很大。第二种方法是通过卷积来求。这种方法直接向我们表明了傅立叶变换的实质。让一个频率为f正弦基信号沿着输入信号进行卷积，检测与其频率相同的子信号。我们知道，在一个周期内，正弦曲线有个性质：同频率的正弦函数相乘积分为1，不同频率的则为0。(上面是个不太准确的定义，严格定义请参考高等数学)。这样，‘游走’一圈之后，对应频率位值的值就会非常大，而其他地方则为0.(这是理想的结果，现实中的信号多带有噪声，所以结果有些出入)。第三种方法就是顶顶大名的快速傅立叶变换了(FFT)。这个方法可以把计算时间复杂度维持在nlog(n)数量级上，这在数据量非常大的时候，速度能够提高成千倍。FFT的思想早在高斯的时候就诞生了，不过直到1967年的时候才由两位牛人利用计算机实现了。如果当时有计算机的话，估计就轮不到这两哥们了。关于FFT的计算方法，我将在陆续给出详细解释(说实话，我到现在还没有完全弄懂，所以就不误导大家了)。
OK，废话说了太多了，现在来分析一下这个小程序吧。这个程序是利用CImg图像库做的。CImg是一个不错的图像库(当然，它也可以用来处理1D信号了)，该库非常的轻巧(就一个头文件，使用的时候include进去就可以了)，该库由于利用的模板编程的方法，所以使得该程序具有良好的可重用性。程序的前几行，主要是定义了一些常量，紧接着就是初始化信号，这个信号具有50Hz的频率和50.0的幅值，信号的采样频率为1024(这个为了FFT， 2的次方)。接下来，就是把我们的输入信号画出来。然后，利用get_FFT函数进行FFT(参数’y‘表明是在y方向进行1D变换)。最后，求取FFT的Magitude。
在结果图像中，0的位置为直流分量(工程师都这么叫)，其实就是输入信号的均值。我们发现在50的位置，值特别大(这里我加了缩放因子的，否则显示出来不太美观)。眼力好的人或许发现这个图像貌似的是对称。恩，是的。所以在显示的时候，我们可以只显示一半。从图中，我们可以知道，输入信号有频率为50的正弦信号组成。实际信号中，可能会有几个峰值，最大的叫first frequency,其次呢?当然是secocnd frequency喽，呵呵。Magitude用来表示能量，那相位(Phase)呢?它是用来表示正弦信号的初始位置的，比如cos(2×PI×f + 45)，那么这个信号的初始相位就为45(这里的单位是度，而函数库中所提够的算法结果一般为弧度)。由于我们的信号初始相位为0，所以你看到第三幅图像中所有的值均为0.
转自：http://blog.sina.com.cn/s/blog_61d0dee10100godv.html
通过FFT计算实例，进一步了解蝶形图
家庭作业式的例题，实际计算FFT中两种算法DIT和DIF，证明两者的结果完全相同。更好地了解蝶形图中的记号和术语。最好将蝶形图和计算步骤式对照来读。如果有笔误也能看出来。其中DIT图已经将计算结果注在图上了。DIF则没有注在图上。转自：http://blog.sina.com.cn/s/blog_3feefc7c0102vjgn.html
谱密度、功率谱密度和能量谱密度
当波的频谱密度乘以一个适当的系数后将得到每单位频率波携带的功率，这被称为信号的功率谱密度（power spectral density, PSD）或者谱功率分布（spectral power distribution, SPD）。功率谱密度的单位通常用每赫兹的瓦特数（W/Hz）表示，或者使用波长而不是频率，即每纳米的瓦特数（W/nm）来表示。能量谱密度能量谱密度描述的是信号或者时间序列（应该就是我们平时所说的随时间而变的信号或者函数或者物理量）的能量或者变化如何随着频率分布。如果 f(t) 是一个有限能量信号，即平方可积，那么信号的谱密度 Φ(ω) 就是信号连续傅里叶变换幅度（体现从时域到频域的一种变化幅度，在时域中变化越快表明周期越短，频率约大，那么变化到频域中也应该有对应的特征）的平方。其中 ω 是角频率（循环频率的 2π 倍），F(ω) 是 f(t) 的连续傅里叶变换。 F *(ω)是F(ω)的共轭函数。如果信号是离散的 fn，经过有限的元素之后，仍然得到能量谱密度：其中 F(ω) 是 fn 的离散时间傅里叶变换。如果所定义的数值个数是有限的，这个序列可以看作是周期性的，使用离散傅里叶变换得到离散频谱，或者用零值进行扩充从而可以作为无限序列的情况计算谱密度。乘数因子 1 / 2π 经常不是绝对的，它随着不同傅里叶变换定义的归一化常数的不同而不同。 功率谱密度上面能量谱密度的定义要求信号的傅里叶变换必须存在，也就是说信号平方可积或者平方可加。一个经常更加有用的替换表示是功率谱密度（PSD），它定义了信号或者时间序列的功率如何随频率分布。这里功率可能是实际物理上的功率，或者更经常便于表示抽象的信号被定义为信号数值的平方，也就是当信号的负载为1欧姆(ohm)时的实际功率。此瞬时功率（平均功率的中间值）可表示为：由于平均值不为零的信号不是平方可积的，所以在这种情况下就没有傅里叶变换。幸运的是维纳-辛钦定理（Wiener-Khinchin theorem）提供了一个简单的替换方法，如果信号可以看作是平稳随机过程，那么功率谱密度就是信号自相关函数的傅里叶变换。信号的功率谱密度当且仅当信号是广义的平稳过程的时候才存在。如果信号不是平稳过程，那么自相关函数一定是两个变量的函数，这样就不存在功率谱密度，但是可以使用类似的技术估计时变谱密度。转自：http://blog.sina.com.cn/s/blog_6419abc70100x8us.html快速傅立叶变换的问题
快速傅氏变换，是离散傅氏变换的快速算法，它是根据离散傅氏变换的奇、偶、虚、实等特性，对离散傅立叶变换的算法进行改进获得的。它对傅氏变换的理论并没有新的发现，但是对于在计算机系统或者说数字系统中应用离散傅立叶变换，可以说是进了一大步。设x(n)为N项的复数序列，由DFT变换，任一X(m)的计算都需要N次复数乘法和N-1次复数加法，而一次复数乘法等于四次实数乘法和两次实数加法，一次复数加法等于两次实数加法，即使把一次复数乘法和一次复数加法定义成一次“运算”(四次实数乘法和四次实数加法)，那么求出N项复数序列的X(m), 即N点DFT变换大约就需要N2次运算。当N=1024点甚至更多的时候，需要N2=1048576次运算，在FFT中，利用WN的周期性和对称性，把一个N项序列(设N=2k,k为正整数)，分为两个N/2项的子序列，每个N/2点DFT变换需要(N/2)2次运算，再用N次运算把两个N/2点的DFT 变换组合成一个N点的DFT变换。这样变换以后，总的运算次数就变成N+2(N/2)2=N+N2/2。继续上面的例子，N=1024时，总的运算次数就变成了525312次，节省了大约50%的运算量。而如果我们将这种“一分为二”的思想不断进行下去，直到分成两两一组的DFT运算单元，那么N点的 DFT变换就只需要Nlog2N次的运算，N在1024点时，运算量仅有10240次，是先前的直接算法的1%，点数越多，运算量的节约就越大，这就是 FFT的优越性.
傅里叶变换(Transformée de Fourier)是一种积分变换。因其基本思想首先由法国学者傅里叶系统地提出，所以以其名字来命名以示纪念。
应用
傅里叶变换在物理学、数论、组合数学、信号处理、概率论、统计学、密码学、声学、光学、海洋学、结构动力学等领域都有着广泛的应用(例如在信号处理中，傅里叶变换的典型用途是将信号分解成幅值分量和频率分量)。
概要介绍
傅里叶变换能将满足一定条件的某个函数表示成三角函数(正弦和/或余弦函数)或者它们的积分的线性组合。在不同的研究领域，傅里叶变换具有多种不同的变体形式，如连续傅里叶变换和离散傅里叶变换。最初傅里叶分析是作为热过程的解析分析的工具被提出的(参见：林家翘、西格尔著《自然科学中确定性问题的应用数学》，科学出版社，北京。原版书名为 C. C. Lin & L. A. Segel, Mathematics Applied to Deterministic Problems in the Natural Sciences, Macmillan Inc., New York, 1974)。
傅里叶变换属于谐波分析。
傅里叶变换的逆变换容易求出,而且形式与正变换非常类似;
正弦基函数是微分运算的本征函数,从而使得线性微分方程的求解可以转化为常系数的代数方程的求解.在线性时不变的物理系统内,频率是个不变的性质,从而系统对于复杂激励的响应可以通过组合其对不同频率正弦信号的响应来获取;
卷积定理指出:傅里叶变换可以化复杂的卷积运算为简单的乘积运算,从而提供了计算卷积的一种简单手段;
离散形式的傅里叶变换可以利用数字计算机快速的算出(其算法称为快速傅里叶变换算法(FFT)).
转自：http://blog.sina.com.cn/s/blog_4b2f8f9601008zqh.html
页:
[1]