fft原理详解Word下载.docx
《fft原理详解Word下载.docx》由会员分享，可在线阅读，更多相关《fft原理详解Word下载.docx（48页珍藏版）》请在冰点文库上搜索。
用数学表达式就是如下：
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的结果的模值如图所示。
图1FFT结果
从图中我们可以看到，在第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。
接着，我们来计算各点的幅度值。
分别计算这三个点的模值，
结果如下：
512
384
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程序]
实例一：
S=2+3cos（2pi*50t-pi/6）+1.5cos（2pi*75t+pi/2）
closeall;
%先关闭所有图片
Adc=2;
%直流分量幅度
A1=3;
%频率F1信号的幅度
A2=1.5;
%频率F2信号的幅度
F1=50;
%信号1频率（Hz）
F2=75;
%信号2频率（Hz）
Fs=256;
%采样频率（Hz）
P1=-30;
%信号1相位（度）
P2=90;
%信号相位（度）
N=256;
%采样点数
t=[0:
1/Fs:
N/Fs];
%采样时刻
%信号
S=Adc+A1*cos（2*pi*F1*t+pi*P1/180）+A2*cos（2*pi*F2*t+pi*P2/180）;
%显示原始信号
plot（S）;
title（'
原始信号'
）;
figure;
Y=fft（S,N）;
%做FFT变换
Ayy=（abs（Y））;
%取模
plot（Ayy（1:
N））;
%显示原始的FFT模值结果
FFT模值'
Ayy=Ayy/（N/2）;
%换算成实际的幅度
Ayy
（1）=Ayy
（1）/2;
F=（[1:
N]-1）*Fs/N;
%换算成实际的频率值
plot（F（1:
N/2）,Ayy（1:
N/2））;
%显示换算后的FFT模值结果
幅度-频率曲线图'
Pyy=[1:
N/2];
fori=1:
N/2
Pyy（i）=phase（Y（i））;
%计算相位
Pyy（i）=Pyy（i）*180/pi;
%换算为角度
end;
N/2）,Pyy（1:
%显示相位图
相位-频率曲线图'
S=1+0.1cos（2pi*20t）+0.2cos（2pi*60t）
Adc=1;
A1=0.1;
A2=0.2;
F1=20;
F2=60;
P1=0;
P2=0;
实例三
Fs=512;
N=1024;
（N-1）/Fs];
%%显示原始信号
%%FFT变换后
%%幅度频率曲线图
%%相位频率曲线图
fori=1:
Pyy（i）=phase（Y（i））;
Pyy（i）=Pyy（i）*180/pi;
实例四
关于FFT的相位谱
（2011-07-1311:
41:
56）
转载▼
标签：
相位谱正弦信号延拓进行it
分类：
机械技术
先看一下我收到的程序，作为研究对象的信号是这样产生的：
T=128;
N=128;
dt=T/N;
t=dt*（1:
N）;
x=2*cos（2*t-pi/4）;
...
（我觉得这个信号存在一点问题，因为t是从1开始的，所以它的初相应该和-pi/4有点差别吧。
）
为什么进行FFT，用angle得到相位－频率特性却不能反映这个信号的初始相位？
胡广书的《数字信号处理－理论、算法与实现（第二版）》第三章第八节《关于正弦信号抽样的讨论》，得出了关于正弦信号抽样的六个结论，最后总结了一个总的原则：
抽样频率应为信号频率的整数倍，抽样点数应包含整周期。
书中的结论五与采样频率和抽样点数有很大的关联。
结论五主要说只有满足了上面的那个总的原则，频谱泄漏才不会发生。
我想不光是幅度谱的频谱泄漏现象，抽样频率和抽样点数同样会对相位谱产生影响。
考虑一个无限长的正弦信号（相当于初相为－90°
），如果我们截取它的整数个周期，然后对截短的信号进行周期延拓，则得到的延拓的信号与原无限长正弦没有区别。
现在我们再次对这个无限长的正弦进行截短，长度为1.5个周期，然后对截短信号进行周期延拓，看看我们得到了什么？
下图、无限长正弦：
下图、截短信号
下图、对截短信号周期延拓：
可以看出，此时进行周期延拓得到的信号与原来的正弦信号大相径庭。
新的周期信号是一个周期的偶函数，原无限长正弦是一个周期的奇函数，两者奇偶性都不一样了，因此不能指望利用新的信号的DFT求出原信号的初相。
exp（-jωt）=cos（ωt）-jsin（ωt），进行变换的时候，若f（t）为实偶函数，则f（t）sin（ωt）就是奇函数，对一个奇函数在对称区间内积分只能得到0，因此实偶函数的傅立叶变换肯定是实的，对一个实数用angle求相位，当然相位是0。
而原正弦肯定是初相为－90°
。
我想这就是问题所在，DFT就是DFS，只不过DFT先将有限长信号进行周期延拓，然后求DFS，再截取一个周期。
使用DFT，在有限的观测时间内采集信号的信息。
如若观测时间内正好得到了整数个正弦周期，则DFT的周期延拓可以不失真的表示原正弦，可是如果观测时间内得到的信号不是整数个周期，那么问题随之而来，就像上面的例子，观测时间内得到了1.5个周期的正弦，然后进行周期延拓，显然乱了套。
如果满足了胡广书老师所总结的抽样条件，则对正弦的DFT谱无疑可以很好地反映初相，我写了两个例子：
第一个例子，信号只包含一个正弦：
t=linspace（0,2-0.125,16）;
x=cos（2*pi*t+pi/4）;
X=fft（x）;
stem（abs（X））;
stem（angle（X）/pi*180）;
幅度谱：
相位谱：
可以看见DFT相位谱第三个点对应正弦的相位，刚好是45°
第二个例子信号中包含两个正弦：
x=cos（2*pi*t+pi/4）+2*cos（2*pi*0.5*t+pi/8）;
可以看见DFT相位谱第二个和第三个点对应两个正弦的相位，刚好是22.5°
和45°
如果没有满足上面所说的条件，就会得到不准确的结果，有兴趣可试试下面的代码：
t=linspace（0,2.5-0.125,32）;
如何克服这个问题？
我觉得这非常困难。
在不能预知信号频率的情况下，无法确定采样频率和观测点数。
也许可以先进行一次观测，通过幅度谱估计出正弦的频率，然后根据频率调整抽样频率，重新对信号进行采样，使采样符合上面所述的条件。
但是这样做有很多的问题，例如硬件可能不好实现。
而且虽然第二次调整了采样频率和抽样点数，可是初始相位已经无法得到了，因为第二次采样不可能再从零时刻开始。
Sandygreta同学说可以这样做，先以较高的抽样频率对信号进行采样，通过FFT幅度谱估算出正弦信号的频率，然后计算出满足抽样条件的最佳的抽样频率和观测时间，使抽样频率为正弦频率的整数倍（大于2倍），且观测时间内能正好得到整数个正弦周期。
然后对刚才采集的信号样本进行插值，接着使用计算出来的采样频率和观测时间对插值的结果重新采样，计算FFT，得到初始相位。
实例五
x=0:
.001:
1;
fs=1000;
y=sin（2*pi*50*x）;
N=64;
fork=1:
4
N=N*2;
M=fft（y,N）;
Py=abs（M）*2/N;
f=fs*（0:
N/2）/N;
%fs采样率
subplot（4,1,k）;
stem（f,Py（1:
N/2+1））;
xlim（[0100]）;
title（['
N='
num2str（N）]）;
end
实例六
同频率信号滤波的问题
发布时间：
2011-06-1417:
36
技术类别：
测试测量
aiyou18：
小子我有一个振动发生器，产生50Hz的振动信号。
这个振动信号作用到一个小球上，让这个小球在一个一个台面上摩擦。
摩擦力由力传感器测出。
理论上这个摩擦力应该是个50HZ的方波信号。
现在的问题是，振动发生器也固定在台面上，因此台面受振动发生器反作用，产生一个50hz的正弦波，这个信号也被叠加到摩擦力信号里。
我想问的是，同样是50HZ的正弦波和50HZ的方波叠加到一起，如何把正弦波隔离掉只留下50hz的方波？
他们之间的相位差不确定。
-------------------------------------
可以根据100hz处的谐波分量推算方波的相位和幅值
------------------------------------
可是,我那个方波不是标准的方波信号，类似于梯形波，也就是波形顶端是个斜面。
如果采用计算机采样，我如何把100hz的谐波取出呢？
-----------------------------------
只要你那个正弦波是标准的正弦波就可以用高次谐波的信息。
找一下梯形波的傅立叶级数展开式，找到高次谐波参数与基波之间的关系（有影响参数肯定有斜面的坡度，占空比等，这些当然你需要事先知道）。
对采集的信号（最好整周期采样），做FFT，在谱图上第二个明显谱峰的地方就应该是二次谐波。
如果你的正弦是标准的，那么这个谱峰就是梯形波参数相关。
根据它能推出参数来。
我不记得有直接公式，需要你自己弄。
弄好了之后，也许写一篇小文章没问题。
--------------------------------------
谢谢Master。
确实同频率的信号滤波，我还第一次遇到。
而且还是真实使用中的压电型力传感器产生的，因此毛刺干扰都很多。
查过很多资料，万方的数据库中，也没有相关的资料供借鉴。
幸好你给我提供了一个思路。
我会去试试，不明白再求教。
再问：
我是以5k采样率对一个波形采集，也就是间隔0.2ms保存一个数值。
最后得到50个数值的这样一组数据，我想用matlab得到这个数组的频率，幅值曲线。
请问如何编程？
data.txt中的一列数据是：
431
827
591
474......
matlab程序我是这么写的：
a=load（'
c:
\data.txt'
y=abs（fft（a））
plot（y）
这个时候出来一个曲线，我不知道x轴表示什么意思?
怎么才能让x轴是hz单位？
直接画图是数据序列号，第一个点对于0频率，以后的增量是1/采样长度。
这是方波谱图
这是锯齿波谱图
真实时域曲线
这个是我采集来的真实曲线。
这个是我对真实时域曲线做的谱图
这个是我们想要的理想曲线（时域）
这个是我们想要的理想曲线（时域
这个是对时域曲线作的频谱图
我们使用的试验机
情况在图上基本已经注明了：
激振器产生一个50hz的正弦波，推动一个金属小球，在一个“派”型支架上摩擦。
派型支架，把摩擦力传递到PCB传感器上。
我们的摩擦力曲线是由PCB力传感器测量出的。
根据我们的理论分析，摩擦力应该是50hz的，有一个小尖的，梯形波。
但是实际过程中，激振器会对底座有反作用力，产生一个50hz的正弦振动。
这个振动也会被PCB传感器测出。
这个正弦波是我们不想要的。
使得实际测量的波形类似一个锯齿波。
对真实波形做FFT，发现频谱图中100HZ的分量很高。
我弄不懂，一个50HZ的振动系统里，怎么会冒出这么多的100hz的波形的？
对方波做FFT，发现方波中根本不含100HZ分量。
都是50，150，250等50的奇数倍频率分量。
而干扰振动信号应该是50hz正弦波。
我想问的是：
1.这100hz信号可能是怎么产生的？
如何滤掉。
2.如何把50hz的底座的正弦振动滤掉一部分？
但是摩擦力波形里50hz分量需要保留。
----------------------------------------
yangzj：
方波是只有1，3，5等奇次谐波，但你的信号是类似于你说的锯齿波呀，他的谐波你也画出来了，各个整数次谐波都有。
vibrationmaster的意思是说，如果知道理想曲线的时域表达式，那就可以用傅立叶级数公式求出它的各个谐次的表达式，这些谐次的幅值和相位都是相互关联的，那么就可以通过50Hz以外的谐波幅值和相位推算出这个信号50Hz处的幅值和相位，这样就可去掉正弦振动的影响。
另外，50Hz与工频重合了，工频及其谐波会有工扰，为什么不把激振频率与工频错开呢？
还有应该做整周期采样结果才会准。
----------------------------
toyangzi：
你没明白我的意思。
那个锯齿波是我实际采集来的信号。
可是这个信号的形状与我预期的不一样。
我预期的信号波形，也就是理论分析得出的波形，应该是那个有一个小尖的，梯形的波形。
我不明白为什么实际波形和理论波形差的那么远？
其中有一部分干扰是50hz的正弦波。
这个是我知道的。
另一部份的100hz分量，我不晓得是哪里来的。
我想得到一个滤波算法，或者别的什么实时处理方法，把我的锯齿波，变成类似于梯形波的形状。
当然这个转变是要有据可循的。
要有理论依据。
我之所以采用50HZ频率的激振器做实验。
是因为人家老外是这么干的，他们立了行业标准！
我没办法改变。
我不知道理想曲线的时域表达式。
升级会员 