利⽤matlab怎样进⾏频谱分析
图像的频率是表征图像中灰度变化剧烈程度的指标，是灰度在平⾯空间上的梯度。如：⼤⾯积的沙漠在图像中是⼀⽚灰度变化缓慢的区域，对应的频率值很低；⽽对于地表属性变换剧烈的边缘区域在图像中是⼀⽚灰度变化剧烈的区域，对应的频率值较⾼。傅⽴叶变换在实际中有⾮常明显的物理意义，设f是⼀个能量有限的模拟信号，则其傅⽴叶变换就表⽰f的谱。从纯粹的数学意义上看，傅⽴叶变换是将⼀个函数转换为⼀系列周期函数来处理的。从物理效果看，傅⽴叶变换是将图像从空间域转换到频率域，其逆变换是将图像从频率域转换到空
间域。换句话说，傅⽴叶变换的物理意义是将图像的灰度分布函数变换为图像的频率分布函数，傅⽴叶逆变换是将图像的频率分布函数变换为灰度分布函数。
这样通过观察傅⽴叶变换后的频谱图，也叫功率图，我们⾸先就可以看出，图像的能量分布，如果频谱图中暗的点数更多，那么实际图像是⽐较柔和的（因为各点与邻域差异都不⼤，梯度相对较⼩），反之，如果频谱图中亮的点数多，那么实际图像⼀定是尖锐的，边界分明且边界两边像素差异较⼤的。对频谱移频到原点以后，可以看出图像的频率分布是以原点为圆⼼，对称分布的。将频谱移频到圆⼼除了可以清晰地看出图像频率分布以外，还有⼀个好处，它可以分离出有周期性规律的⼲扰信号，⽐如正弦⼲扰，⼀副带有正弦⼲扰，移频到原点的频谱图上可以看出除了中⼼以外还存在以某⼀点为中⼼，对称
分布的亮点集合，这个集合就是⼲扰噪⾳产⽣的，这时可以很直观的通过在该位置放置带阻滤波器消除⼲扰。另外我还想说明以下⼏点：
1、图像经过⼆维傅⽴叶变换后，其变换系数矩阵表明：
若变换矩阵Fn原点设在中⼼，其频谱能量集中分布在变换系数短阵的中⼼附近(图中阴影区)。若所⽤的⼆维傅⽴叶变换矩阵Fn的原点设在左上⾓，那么图像信号能量将集中在系数矩阵的四个⾓上。这是由⼆维傅⽴叶变换本⾝性质决定的。同时也表明⼀股图像能量集中低频区域。
2、变换之后的图像在原点平移之前四⾓是低频，最亮，平移之后中间部分是低频，最亮，亮度⼤说明低频的能量⼤（幅⾓⽐较⼤）。
从计算机处理精度上就不难理解，⼀个长度为N的信号，最多只能有N/2+1个不同频率，再多的频率就超过了计算机所能所处理的精度范围）
X[]数组⼜分两种，⼀种是表⽰余弦波的不同频率幅度值：Re X[]，另⼀种是表⽰正弦波的不同频率幅度值：Im X[]，Re是实数(Real)的意思，Im是虚数(Imagine)的意思，采⽤复数的表⽰⽅法把正余弦波组合起来进⾏表⽰，但这⾥我们不考虑复数的其它作⽤，只记住是⼀种组合⽅法⽽已，⽬的是为了便于表达（在后⾯我们会知道，复数形式的傅⽴叶变换长度是N，⽽不是N/2+1）。
⽤Matlab实现快速傅⽴叶变换
FFT是离散傅⽴叶变换的快速算法，可以将⼀个信号变换到频域。有些信号在时域上是很难看出什么特征的，但是如果变换到频域之后，就很容易看出特征了。这就是很多信号分析采⽤FFT变换的原因。另外，FFT可以将⼀个信号的频谱提取出来，这在频谱分析⽅⾯也是经常⽤的。
虽然很多⼈都知道FFT是什么，可以⽤来做什么，怎么去做，但是却不知道FFT之后的结果是什意思、如何决定要使⽤多少点来做FFT。
现在就根据实际经验来说说FFT结果的具体物理意义。⼀个模拟信号，经过ADC采样之后，就变成了数字信号。采样定理告诉我们，采样频率要⼤于信号频率的两倍，这些我就不在此啰嗦了。
采样得到的数字信号，就可以做FFT变换了。N个采样点，经过FFT之后，就可以得到N个点的FFT结果。为了⽅便进⾏FFT运算，通
常N取2的整数次⽅。
假设采样频率为Fs，信号频率F，采样点数为N。那么FFT之后结果就是⼀个为N点的复数。每⼀个点就对应着⼀个频率点。这个点的模值，就是该频率值下的幅度特性。具体跟原始信号的幅度有什么关系呢？假设原始信号的峰值为A，那么FFT的结果的每个点（除了第⼀个点直流分量之外）的模值就是
A的N/2倍。⽽第⼀个点就是直流分量，它的模值就是直流分量的N倍。⽽每个点的相位呢，就是在该频率下的信号的相位。第⼀个点表⽰直流分量（即0Hz），⽽最后⼀个点N的再下⼀个点（实际上这个点是不存在的，这⾥是假设的第N+1个点，也可以看做是将第⼀个点分做两半分，另⼀半移到最后）则表⽰采样频率Fs，这中间被N-1个点平均分成N等份，每个点的频率依次增加。例如某点n所表⽰的频率为：Fn=(n-1)*Fs/N。由上⾯的公式可以看出，Fn所能分辨到频率为Fs/N，如果采样频率Fs为1024Hz，采样点数
为1024点，则可以分辨到1Hz。1024Hz的采样率采样1024点，刚好是1秒，也就是说，采样1秒时间的信号并做FFT，则结果可以分析
到1Hz，如果采样2秒时间的信号并做FFT，则结果可以分析到0.5Hz。如果要提⾼频率分辨⼒，则必须增加采样点数，也即采样时间。频率分辨率和采样时间是倒数关系。
假设FFT之后某点n⽤复数a+bi表⽰，那么这个复数的模就是An=根号a*a+b*b，相位就是Pn=atan2(b,a)。根据以上的结果，就可以计算
出n点（n≠1，且n<=N/2）对应的信号的表达式为：An/(N/2)*cos(2*pi*Fn*t+Pn)，即2*An/N*cos(2*pi*Fn*t+Pn)。对于n=1点的信号，是直流分量，幅度即为A1/N。由于FFT结果的对称性，通常我们只使⽤前半部分的结果，即⼩于采样频率⼀半的结果。
下⾯以⼀个实际的信号来做说明。假设我们有⼀个信号，它含有2V的直流分量，频率为50Hz、相位为-30度、幅度为3V的交流信号，以及⼀个频率（f0）为75Hz、相位为90度、幅度为1.5V的交流信号。⽤数学表达式就是如下：S=2+3*cos(2*pi*50*t-
pi*30/180)+1.5*cos(2*pi*75*t+pi*90/180)。式中cos参数为弧度，所以-30度和90度要分别换算成弧度。我们以256Hz的采样率对这个信号进⾏采样，总共采样256点。按照我们上⾯的分析，Fn=(n-1)*Fs/N，我们可以知道，每两个点之间的间距就是1Hz，第n个点的频率就是n-1。我们的信号有3个频率：0Hz、50Hz、75Hz，应该分别在第1个点、第51个点、第76个点上出现峰值，其它各点应该接近0。实际情况如何呢？我们来看看FFT的结果的模值如图所⽰。
从图中我们可以看到，在第1点、第51点、和第76点附近有⽐较⼤的值。我们分别将这三个点附近的数据拿上来细看：
1点：512+0i
2点：-2.6195E-14 - 1.4162E-13i
3点：-2.8586E-14 - 1.1898E-13i
50点：-6.2076E-13 - 2.1713E-12i
51点：332.55 - 192i
52点：-1.6707E-12 - 1.5241E-12i
75点：-2.2199E-13 -1.0076E-12i
76点：3.4315E-12 + 192i
77点：-3.0263E-14 +7.5609E-13i
很明显，1点、51点、76点的值都⽐较⼤，它附近的点值都很⼩，可以认为是0，即在那些频率点上的信号幅度为0。接着，我们来计算各点的幅度值。分别计算这三个点的模值，结果如下：
1点：512
51点：384
76点：192
按照公式，可以计算出直流分量为：512/N=512/256=2；50Hz信号的幅度为：384/(N/2)=384/(256/2)=3；75Hz信号的幅度
为192/(N/2)=192/(256/2)=1.5。可见，从频谱分析出来的幅度是正确的。
然后再来计算相位信息。直流信号没有相位可⾔，不⽤管它。先计算50Hz信号的相位，atan2(-192, 332.55)=-0.5236，结果是弧度，换算为⾓度就是180*(-0.5236)/pi=-30.0001。再计算75Hz信号的相位，atan2(192, 3.4315E-12)=1.5708弧度，换算成⾓度就
是180*1.5708/pi=90.0002。可见，相位也是对的。根据FFT结果以及上⾯的分析计算，我们就可以写出信号的表达式了，它就是我们开始提供的信号。
总结：假设采样频率为Fs，采样点数为N，做FFT之后，某⼀点n（n从1开始）表⽰的频率为：Fn=(n-1)*Fs/N；该点的模值除以N/2就是对应该频率下的信号的幅度（对于直流信号是除以N）；该点的相位即是对应该频率下的信号的相位。相位的计算可⽤函数atan2(b,a)计算。atan2(b,a)是求坐标为(a,b)点的⾓度值，范围从-pi到pi。要精确到xHz，则需要采样长度为1/x秒的信号，并做FFT。要提⾼频率分辨率，就需要增加采样点数，这在⼀些实际的应⽤中是不现实的，需要在较短的时间内完成分析。解决这个问题的⽅法有频率细分法，⽐较简单的⽅法是采样⽐较短时间的信号，然后在后⾯补充⼀定数量的0，使其长度达到需要的点数，再做FFT，这在⼀定程度上能够提⾼频率分辨⼒。具体的频率细分法可参考相关⽂献。
附贴上上述例⼦的matlab程序：
Matlab的例⼦（⼀）
t=0:1/256:1;%采样步长
y= 2+3*cos(2*pi*50*t-pi*30/180)+1.5*cos(2*pi*75*t+pi*90/180);
N=length(t); %样点个数
plot(t,y);
fs=256;%采样频率
df=fs/(N-1);%分辨率
f=(0:N-1)*df;%其中每点的频率
Y=fft(y)/N*2;%真实的幅值
%Y=fftshift(Y);
figure(2)
plot(f,abs(Y));
由于以上程序是结合傅⾥叶算法转换得到的对称图，⽽常⽤的只需要⼀半就可以了。对应的程序如下：
t=0:1/256:1;%采样步长
y= 2+3*cos(2*pi*50*t-pi*30/180)+1.5*cos(2*pi*75*t+pi*90/180);
N=length(t); %样点个数
plot(t,y);
fs=256;%采样频率
df=fs/(N-1);%分辨率
f=(0:N-1)*df;%其中每点的频率matlab软件怎么使用
Y=fft(y(1:N))/N*2;%真实的幅值
%Y=fftshift(Y);
figure(2)
plot(f(1:N/2),abs(Y(1:N/2)));
由于⽔平有限，难免会出现错误之处，欢迎指正，恳请批评！

版权声明：本站内容均来自互联网，仅供演示用，请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系QQ:729038198，我们将在24小时内删除。