fft-原理详解.doc

FFT算法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点（n1，且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的交流信号，以及一个频率为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 FFT结果 从图中我们可以看到，在第1点、第51点、和第76点附近有比较大的值。我们分别将这三个点附近的数据拿上来细看：1点： 512+0i2点： -2.6195E-14 - 1.4162E-13i3点： -2.8586E-14 - 1.1898E-13i50点：-6.2076E-13 - 2.1713E-12i51点：332.55 - 192i52点：-1.6707E-12 - 1.5241E-12i75点：-2.2199E-13 -1.0076E-12i76点：3.4315E-12 + 192i77点：-3.0263E-14 +7.5609E-13i 很明显，1点、51点、76点的值都比较大，它附近的点值都很小，可以认为是0，即在那些频率点上的信号幅度为0。接着，我们来计算各点的幅度值。分别计算这三个点的模值，结果如下：1点： 51251点：38476点：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)close all; %先关闭所有图片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模值结果title(FFT 模值);figure;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模值结果title(幅度-频率曲线图);figure;Pyy=1:N/2;for i=1:N/2Pyy(i)=phase(Y(i); %计算相位Pyy(i)=Pyy(i)*180/pi; %换算为角度end;plot(F(1:N/2),Pyy(1:N/2); %显示相位图title(相位-频率曲线图);实例一：S=1+0.1cos(2pi*20t)+0.2cos(2pi*60t)close all; %先关闭所有图片Adc=1; %直流分量幅度A1=0.1; %频率F1信号的幅度A2=0.2; %频率F2信号的幅度F1=20; %信号1频率(Hz)F2=60; %信号2频率(Hz)Fs=256; %采样频率(Hz)P1=0; %信号1相位(度)P2=0; %信号相位(度)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模值结果title(FFT 模值);figure;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模值结果title(幅度-频率曲线图);figure;Pyy=1:N/2;for i=1:N/2Pyy(i)=phase(Y(i); %计算相位Pyy(i)=Pyy(i)*180/pi; %换算为角度end;plot(F(1:N/2),Pyy(1:N/2); %显示相位图title(相位-频率曲线图);实例三close all; %先关闭所有图片Adc=2; %直流分量幅度A1=3; %频率F1信号的幅度F1=50; %信号1频率(Hz)P1=-30; %信号1相位(度)A2=1.5; %频率F2信号的幅度F2=75; %信号2频率(Hz)P2=90; %信号相位(度)Fs=512; %采样频率(Hz)N=1024; %采样点数t=0:1/Fs:(N-1)/Fs; %采样时刻%信号S=Adc+A1*cos(2*pi*F1*t+pi*P1/180)+A2*cos(2*pi*F2*t+pi*P2/180);% 显示原始信号plot(S);title(原始信号);% FFT变换后figure;Y = fft(S,N); %做FFT变换Ayy = (abs(Y); %取模plot(Ayy(1:N); %显示原始的FFT模值结果title(FFT 模值);% 幅度频率曲线图figure;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模值结果title(幅度-频率曲线图);% 相位频率曲线图figure;Pyy=1:N/2;for i=1:N/2 Pyy(i)=phase(Y(i); %计算相位 Pyy(i)=Pyy(i)*180/pi; %换算为角度end;plot(F(1:N/2),Pyy(1:N/2); %显示相位图title(相位-频率曲线图);实例四关于FFT的相位谱 (2011-07-13 11: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(-jt)=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);figure;stem(angle(X)/pi*180);幅度谱：相位谱：可以看见DFT相位谱第三个点对应正弦的相位，刚好是45。第二个例子信号中包含两个正弦：t=linspace(0,2-0.125,16);x=cos(2*pi*t+pi/4)+2*cos(2*pi*0.5*t+pi/8);X=fft(x);stem(abs(X);figure;stem(angle(X)/pi*180);幅度谱：相位谱：可以看见DFT相位谱第二个和第三个点对应两个正弦的相位，刚好是22.5和45。如果没有满足上面所说的条件，就会得到不准确的结果，有兴趣可试试下面的代码：t=linspace(0,2.5-0.125,32);x=cos(2*pi*t+pi/4);X=fft(x);stem(abs(X);figure;stem(angle(X)/pi*180);如何克服这个问题？我觉得这非常困难。在不能预知信号频率的情况下，无法确定采样频率和观测点数。也许可以先进行一次观测，通过幅度谱估计出正弦的频率，然后根据频率调整抽样频率，重新对信号进行采样，使采样符合上面所述的条件。但是这样做有很多的问题，例如硬件可能不好实现。而且虽然第二次调整了采样频率和抽样点数，可是初始相位已经无法得到了，因为第二次采样不可能再从零时刻开始。Sandygreta同学说可以这样做，先以较高的抽样频率对信号进行采样，通过FFT幅度谱估算出正弦信号的频率，然后计算出满足抽样条件的最佳的抽样频率和观测时间，使抽样频率为正弦频率的整数倍（大于2倍），且观测时间内能正好得到整数个正弦周期。然后对刚才采集的信号样本进行插值，接着使用计算出来的采样频率和观测时间对插值的结果重新采样，计算FFT，得到初始相位。实例五x=0:.001:1; fs=1000;y=sin(2*pi*50*x);N=64;for k=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(0 100); title(N= num2str(N);end实例六同频率信号滤波的问题发布时间：2011-06-14 17:41:36技术类别：测试测量aiyou18：小子我有一个振动发生器，产生50Hz的振动信号。这个振动信号作用到一个小球上，让这个小球在一个一个台面上摩擦。摩擦力由力传感器测出。理论上这个摩擦力应该是个50HZ的方波信号。现在的问题是，振动发生器也固定在台面上，因此台面受振动发生器反作用，产生一个50hz的正弦波，这个信号也被叠加到摩擦力信号里。我想问的是，同样是50HZ的正弦波和50HZ的方波叠加到一起，如何把正弦波隔离掉只留下50hz的方波？他们之间的相位差不确定。-可以根据100hz处的谐波分量推算方波的相位和幅值-可是,我那个方波不是标准的方波信号，类似于梯形波，也就是波形顶端是个斜面。如果采用计算机采样，我如何把100hz的谐波取出呢？-只要你那个正弦波是标准的正弦波就可以用高次谐波的信息。找一下梯形波的傅立叶级数展开式，找到高次谐波参数与基波之间的关系（有影响参数肯定有斜面的坡度，占空比等，这些当然你需要事先知道）。对采集的信号（最好整周期采样），做FFT，在谱图上第二个明显谱峰的地方就应该是二次谐波。 如果你的正弦是标准的，那么这个谱峰就是梯形波参数相关。 根据它能推出参数来。我不记得有直接公式，需要你自己弄。 弄好了之后，也许写一篇小文章没问题。-谢谢Master。确实同频率的信号滤波，我还第一次遇到。而且还是真实使用中的压电型力传感器产生的，因此毛刺干扰都很多。查过很多资料，万方的数据库中，也没有相关的资料供借鉴。幸好你给我提供了一个思路。我会去试试，不明白再求教。-再问：我是以5k采样率对一个波形采集，也就是间隔0.2ms保存一个数值。最后得到50个数值的这样一组数据，我想用matlab得到这个数组的频率，幅值曲线。请问如何编程？data.txt中的一列数据是：431827591474 .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与工频重合了，工频及其谐波会有工扰，为什么不把激振频率与工频错开呢？还有应该做整周期采样结果才会准。-to yangzi：你没明白我的意思。那个锯齿波是我实际采集来的信号。可是这个信号的形状与我预期的不一样。我预期的信号波形，也就是理论分析得出的波形，应该是那个有一个小尖的，梯形的波形。我不明白为什么实际波形和理论波形差的那么远？其中有一部分干扰是50hz的正弦波。这个是我知道的。另一部份的100hz分量，我不晓得是哪里来的。我想得到一个滤波算法，或者别的什么实时处理方法，把我的锯齿波，变成类似于梯形波的形状。当然这个转变是要有据可循的。要有理论依据。我之所以采用50HZ频率的激振器做实验。是因为人家老外是这么干的，他们立了行业标准！我没办法改变。我不知道理想曲线的时域表达式。我只知道50hz正弦干扰振动的时域表达式。我能不能这样：1.把实际信号做FFT变换，得到一系列数组A。（在matlab里就是一系列数组）2.把50hz正弦波做FFT变换，得到一系列数组B。3.c=A-B4.把c做反FFT变换这个时候反变换 得到的曲线是不是去掉了50hz干扰后的信号？另外，你说的“应该做整周期采样结果才会准。”是不是说我做FFT的时候，数据应该是周期的整数倍个啊？-yangzj：你所说的梯形波应该也是各个整数倍频都有的，具体你可用傅立叶级数的公式算一算。并不一定要知道理想曲线的确切时域表达式，有个带参数的模型或许就能找出它各个谐次之间的关系。你的方法原理上行的通，而且没必要做FFT和IFFT，直接时域减就行，但是要保证两者同步。这种信号可能根据信号的特征从时域上进行处理更好些。 引用:原帖由aiyou18于 2007-10-22 14:46 发表另外，你说的“应该做整周期采样结果才会准。”是不是说我做FFT的时候，数据应该是周期的整数倍个啊？是的-to yangzi我因为不知道两个波形之间的相位差的问题。所以我才没有在时域里直接相减。想问一下，在频域里相减需要考虑相位差吗？-频域里复数相减和直接在时域里相减是完全一样的。傅立叶变换具有线性性-请问yangzj：数据是周期的整数倍，我现在有点晕，周期指的是哪个信号周期呢比如下面这个例子，figure(1)subplot(211)f=1000;t=0:1/f:5; x=2*sin(6*pi*t)+5*cos(40*pi*t); plot(t,x)title(两个频率成分3Hz和20Hz复合的正弦信号)xlabel(t/s) ;grid onsubplot(212)M=nextpow2(length(x);N=2M;y=fft(x,N);n=0:N-1;q=f*n/N;stem(q,abs(y)/(length(x)*2)title(正弦信号的LXHFFT f=1000)xlabel(频率Hz) ;grid on;我的f应该设成多大呢？要采样多少个数据，结果才会准确？多谢了-yangzj：让两个成分的频率都是频率分辨率f/N的整数倍就行了-xray：一点个人看法看起来这的确是一个有趣的问题，由于表述和理解上可能存在一些误差，我根据自己的理解，把这个问题整理并重新描述一下，如果有错误，还请 aiyou18 和 yangzj 指正已知条件：1. 方波谱图.jpg : 仿真得到的方波频谱图（理想方波频谱）2. 锯齿波频谱.jpg : 仿真得到的锯齿波频谱图（理想锯齿波频谱）3. 原始图形.jpg : 真实时域曲线（实测）4. 真实波形谱图.jpg : 真实频域曲线（根据实测数据计算）5. 理想曲线.jpg : 仿真得到的时域曲线（理想时域曲线）6. 理想模拟波形.jpg : 仿真得到的频域曲线（根据仿真数据计算）关键问题：1. 对实测数据建模（即产生符合 原始图形.jpg 所示的波形）2. 对实测数据进行处理，滤除干扰成分（即产生符合 理想曲线.jpg 所示的波形）对于问题1，aiyou18从物理装置上分析，认为实测信号应该是一个理想梯形波叠加一个正弦波但是，从 真实波形谱图.jpg 来看，这个模型是不成立的，其理由就是频谱图中多出来了一个100Hz左右的正弦波。因此，我猜想，实测信号应该是一个理想梯形波叠加一个锯齿波（50Hz），这样频域图就可以得到合理的解释了（这一点yangzj也指出来了）。对于问题2，如果把信号模型认定为 理想梯形波叠加一个锯齿波，仅仅从频域滤波来说，还是很难分离的，因为二者的频率出现了比较多的混叠，除非可以知道更多的锯齿波信息。解决思路：（仅供参考）方案一：把PCB动态传感器从底座上分离出来，固定在其他支撑点上，或者在PCB动态传感器和底座之间添加防震材料。评价：物理解决方案，一劳永逸。方案二：增加一个实验，不放置金属小球，让激励器空载运行，或者断开PCB动态传感器和支架的连接装置，直接测量激励器的反作用力（即测量锯齿波），从得到更加详细的锯齿波参数。评价：实现简单，但是空载和加载情况下，锯齿波参数可能会略有差别，而且即时测量得到了详细的锯齿波参数，仍然需要后期的信号处理。To yangzj：这个梯形波的确没有偶次谐波分量，从 理想曲线.jpg上看，去掉直流分量以后，图形应该是关于原点对称的，它是一个奇函数，因此不存在偶次谐波分量。-aiyou18：to xray：非常感谢你对我的帖子的关注！喜欢你搞学问的风格，很多大师都是因为兴趣而科研。你的概括很到位，而且指出的关键点一针见血。 方案一：是最好的办法。这点我也很赞同。可是因为一个实验机，它是需要移动的。所以PCB传感器没法固定在别的上面。只能固定在底座上。 根据牛顿第三定律，激振器输出多少力，就得到多少反作用力。因为振动的频率50hz和振幅1mm是确定的，所以小球的运动速度的平均值v球是一定的。 根据动量守恒：m球 *v球=m底座*v底座。我要减小底座的振动，也就是减小振幅，也就是减小v底座。综上所述，我能做的只能是减小m球，加大m底座。现在我的底座已经用很重的实铁做的，不能再重了。方案二：这个不加小球的实验我做了。测出的空载波形是一个很标准的正弦波，而不是你说的锯齿波。这也是我奇怪的地方：一个正弦波+理论上的梯形波我实际得到的波形。而这个波形据你们分析是含锯齿波的。 我以为用真实波形减去，空载时候测出的正弦波就能得到真实的摩擦力曲线。问题是我不知道真实波形和干扰正弦波的相位差。所以我不知道如何相减？！ 按道理说，这个实验机是个线性系统。应该检测到的任何振动都是50HZ的振动信号（因为启振源是50hz）。可是现在出现了高频信号，而且也可以排除它是电磁干扰。所以我估计摩擦力能使得系统变成非线性系统。所以有可能摩擦力能产生锯齿波？！ 我查了相关论文，有些论文略微提到过。但是试验条件，和我的实验机不同。没法直接比较。而对摩擦力在我的实验机中产生的信号，是不是就一定是个梯形波，我可能还得进一步调查。为什么我说它是梯形波，是由经典力学得出的结论。 滑动摩擦力与速度无关；最大静摩擦略大于动摩擦力。或许在高频振动情况下，摩擦力特性会变化？谁知道呢。但是我知道的是这个真实信号里，一定含有正弦波分量。 我想做的就是从实际测量到的信号里，把这个正弦分量（底座振动引起的）去掉！现在我不晓得相位差，或许能从力学分析上得出，两者的相位差为0？ 楼主我觉得你写整段话的时候，注意一下可读性啊看着太累 -zhangnan3509-to zhangnan3509： 真不好意思。写的兴起，忘了断句了。我简单的说吧：1.这个测量摩擦力的实验机，我发了图，在第一页。激振器引起底座的振动，这个振动信号是一个50hz正弦波。这就是xray所说的不放小球时测出的空载信号。2.我测量出的真实信号，混杂了摩擦力信号A与空载的正弦信号B。两者是同频的50赫兹，相位差我不知道。或许可以从运动力学推出相位差，为0或180度这样的特例。3.根据滑动摩擦力与速度无关，以及最大静摩擦力略大于滑动摩擦力，我们可以推导出：理论上的摩擦力波形是：有个小尖的方波。考虑到有些衰减，所以估计是：有个小尖的梯形波。4.第3点的推导，是基于最简单的力学得到的结果。现代摩擦学，可能复杂的多，因此真实摩擦力曲线也有可能不是梯形波。5.目前，我想做的就是：把真实波形 减去 那个正弦波。难处在于：我不知道相位差。6.当然，第5点的简单相减的前提假设是：空载时的正弦波；在把小球放上去，也就是加载后，正弦波的形状不发生变化。具体变不变，还得我这次定的另一个PCB传感器到货后，我测试了才知道。-yangzj：回复 xray 的帖子To yangzj：这个梯形波的确没有偶次谐波分量，从 理想曲线.jpg上看，去掉直流分量以后，图形应该是关于原点对称的，它是一个奇函数，因此不存在偶次谐波分量。对称性只与实部和虚部是不是为0有关吧，和有没有奇偶谐次有关吗？-xray：回复 yangzj 的帖子不好意思，刚才翻了一下书，是我记错了-楼主能不能把你的测试数据放上来-to yangzi：这个是我采集来的真实信号。也就是几个周期，你够用了吗？采样频率是5K。基本上是100个点一个周期。附件data.txt(3.24 KB)补充的是，我只采集了正半轴的波形。负半轴是对称的。-如果是这样的话，可以解释为什么正弦激励会产生100Hz的分量了，主要是因为采集的数据中把正弦信号负半轴的部分给截掉了仿真程序如下：clearfs = 5000;N = 1000;a1 = ones(N/(50*2),50);a2 = repmat(1:50)/50, N/(50*2),1);a3 = repmat(sin(2*pi*50/fs*(1:50), N/(50*2),1);z = zeros(N/(50*2),50);a = a1 z .;b1 = a(:);a = a2 z .;b2 = a(:);a = a3 z .;b3 = a(:);bfft1 = abs(fft(b1-mean(b1);bfft2 = abs(fft(b2-mean(b2);bfft3 = abs(fft(b3-mean(b3);figure(1); hold on;plot(b1, r);plot(b2, g);plot(b3, b);figure(2); hold on;plot(0:fs/N:fs/2-fs/N, bfft1(1:N/2), r)plot(0:fs/N:fs/2-fs/N, bfft2(1:N/2), g)plot(0:fs/N:fs/2-fs/N, bfft3(1:N/2),