原创 同频率信号滤波的问题

2011-6-14 17:40 2743 12 12 分类: 测试测量
 
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/采样长度。
--------------------------------------
这是方波谱图

59843_1247574802k00C.jpg

这是锯齿波谱图

59843_12475748056tzn.jpg

真实时域曲线
这个是我采集来的真实曲线。
59843_124757481209Zw.jpg
 
这个是我对真实时域曲线做的谱图
这个是我对真实时域曲线做的谱图
59843_1247574817G4Re.jpg
 
这个是我们想要的理想曲线(时域)
这个是我们想要的理想曲线(时域
59843_1247574825e1NT.jpg
 
这个是对时域曲线作的频谱图
这个是对时域曲线作的频谱图
59843_1247574828OPYp.jpg
 
我们使用的试验机
情况在图上基本已经注明了:激振器产生一个50hz的正弦波,推动一个金属小球,在一个“派”型支架上摩擦。派型支架,把摩擦力传递到PCB传感器上。
我们的摩擦力曲线是由PCB力传感器测量出的。根据我们的理论分析,摩擦力应该是50hz的,有一个小尖的,梯形波。
但是实际过程中,激振器会对底座有反作用力,产生一个50hz的正弦振动。这个振动也会被PCB传感器测出。这个正弦波是我们不想要的。
使得实际测量的波形类似一个锯齿波。

对真实波形做FFT,发现频谱图中100HZ的分量很高。我弄不懂,一个50HZ的振动系统里,怎么会冒出这么多的100hz的波形的?
对方波做FFT,发现方波中根本不含100HZ分量。都是50,150,250等50的奇数倍频率分量。
而干扰振动信号应该是50hz正弦波。

我想问的是:
1.这100hz信号可能是怎么产生的?如何滤掉。
2.如何把50hz的底座的正弦振动滤掉一部分?但是摩擦力波形里50hz分量需要保留。
59843_1247574833EF51.jpg
 
----------------------------------------
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-B
4.把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 on
subplot(212)
M=nextpow2(length(x));
N=2^M;
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的分量了,主要是因为采集的数据中把正弦信号负半轴的部分给截掉了
仿真程序如下:

clear

fs = 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), 'b')

To aiyou18:
最好把你测量的空载时候的正弦信号数据也放上来。

-------------------------------------
yangzj:
不知道为什么不测负值时的数据,这样前面提到的通过高次谐波恢复基频的方法也难以实现了。
--------------------------------------
xray:
那也不一定,根据我上面的仿真程序,方波只有奇次谐波,而截断的正弦波是50Hz+偶次谐波
-------------------------------------
yangzj:
也不是说不能,只是会使问题复杂化。
-------------------------------------
VibrationMaster:
to:aiyou18
1.你的方波并非标准方波.
2.这个方波形状是否是固定的? 也就是波宽,占空比,那个左上方的那个三角形等? 我的问题最终是是否只要估计方波高度之后,其他一切都确定了?
3. 如果2的问题答案是yes,用相关法也可以.
---------------------------------------
aiyou18:
to yangzi:因为小球是在“派”型支架上来回摩擦。因此对力传感器的作用是,又拉又压的,因此传感器出来的信号是,有正半轴的信号,也有负半轴的信号。两边的信号是对称的,所以我只取了一个半轴的信号。
 
to yangzi:你说的同步平均是什么意思?把各个周期对应的值相累加后,求平均吗?那只能排出随机噪声,没法消除底座的正弦波干扰啊
 
to xray:我认为产生100hz谐波和只采集了半波信号无关。因为采集卡采集量程可以设置成 -5V到+5V或者是0~10V。为了提高分辨率,在正负信号对称的情况下我只取了正半轴信号,这应该是不影响的。
如果你坚持认为有影响。你可以把正半轴的波形,做一个镜像,对称到负半轴。然后做FFT变换,100hz仍然存在。
----------------------------------------
yangzj:
     引用:
原帖由 aiyou18 于 2007-10-29 13:51 发表
to yangzi:你说的同步平均是什么意思?把各个周期对应的值相累加后,求平均吗?那只能排出随机噪声,没法消除底座的正弦波干扰啊
没法消除正弦波,但减小随机噪声对提取信号很有帮助呀。
 
     引用:
原帖由 aiyou18 于 2007-10-29 14:00 发表
to xray:我认为产生100hz谐波和只采集了半波信号无关。因为采集卡采集量程可以设置成 -5V到+5V或者是0~10V。为了提高分辨率,在正负信号对称的情况下我只取了正半轴信号,这应该是不影响的。
如果你坚持认为有 ...
做镜像对信号的连续性会有影响。我的意思是:如果全部采了,那么照你的理论的话,余弦信号只在50Hz对梯形波有影响,但是负半轴被削掉后,削波后的余弦信号会产生偶次谐波,各梯形波的偶次谐波混在一起。
 
     引用:
原帖由 VibrationMaster 于 2007-10-25 12:26 发表
to:aiyou18
1.你的方波并非标准方波.
2.这个方波形状是否是固定的? 也就是波宽,占空比,那个左上方的那个三角形等? 我的问题最终是是否只要估计方波高度之后,其他一切都确定了?
3. 如果2的问题答案是yes,用相 ...
我觉得VibrationMaster的第二个问题很关键,就是你的理论梯形波能否用参数来表达。
-----------------------------------------
xray:
引用:
你可以把正半轴的波形,做一个镜像,对称到负半轴。然后做FFT变换,100hz仍然存在。
这个方法你试过吗?我试了一下,如果我的程序没有错误的话,结果表明100Hz的成分的确消失了,代码如下:

clear
fs = 5000;
load('data.mat')
Nd = length(data);
Newdata = -data(55:end);
data(5:4+length(Newdata)) = data(5:4+length(Newdata)) + Newdata;
datafft = abs(fft(data-mean(data)));
figure(1);
plot(data)
figure(2);
plot((0:fs:fs/2*Nd-fs)/Nd, datafft(1:floor(Nd/2)))

说明:程序采用的数据是前面帖子中上传的数据,我把它保存为mat文件了。
--------------------------------------
VibrationMaster:
虽然存在,但是肯定小多了。半波整流是非线性操作,肯定会产生高次谐波。你最好保留全波
------------------------------------------
aiyou18:
to VibrationMaster: 刚看了你关于一个矢量,标量的帖子。你说无论矢量,标量都只是表达工具,但是我觉得是并非如此。

我那个方波和小尖(也就是你说的小三角)是不确定的。但是占空比,应该是固定的=50%。幅度也不确定。

因为那是我实际做实验采集来的波形。我想通过不同波形的幅度来区分开不同试验品的润滑能力并把测出的电压值转成真实摩擦力大小。
-------------------------------------------
VibrationMaster:
 
1。一般分析方法很难恢复“小尖”的形状,通过一般的频谱分析也很难保持波形拐角地方的准确形状。 如果你对这些地方感兴趣,那么需要对目标图形参数化描述,然后用实际测量的数据优化这些参数的取值。
2。如果研究数学的人可能认为矢量和标量的差异比较大。但是从物理角度这只是描述工具而已,这是我的观点。比如一个人的体型应该用标量描述,还是矢量描述呢? 我想身高,体重,三围这5个参数大概差不多了,这应该是5个标量吧,可是你非要把这五个标量放在一个{},那么他也可以称为矢量。
3。 如果在一组人群中,体重,身高,胸围,腰围都固定了,变化的只有臀围,那时我们最好用一个标量来表示就可以了,这时用矢量很麻烦。
---------------------------------
aiyou18:
to xray:我把负半轴的波形补全了。然后做fft变换还真是没有100hz的分量了,真神奇!都是50hz的奇数倍频率存在。如果是标准的锯齿波如果正负半轴都有,估计也不会有100hz分量。这个奇数倍频率和偶数倍频率是否出现有什么规律吗?





这是我的程序:
fs=5000;
t=0:1/fs:0.16;
a=load('c:\dataDouble.txt');
number=512;
y=fft(a,number);
n=0:length(y)-1;
f=fs*n/length(y);
plot(f,abs(y))


上传的是我补齐负半轴后的波形数据。
附件
 dataDouble.txt (4.18 KB)

 

我刚才又做了试验。如果把锯齿波正负半轴补全,频谱图里也没有100hz的分量。把理想的模拟波形,补全后,也没有100hz分量。

我可以得出结论,100hz纯粹是因为只采集了正半轴造成的。
-----------------------------------------
VibrationMaster:
那就用3次谐波来恢复。奇周期函数的傅立叶系数中的不含二次谐波是很正常的。
-----------------------------------------------
aiyou18:
to vibrationMaster: 如何用3次谐波(150hz)的幅值推测50hz分量的幅值? 有什么对应关系吗  小弟不是干这个的,不太懂。请指教。

另外,请教一个问题,波形A上叠加50hz的正弦波B生成了C。是不是不论A与B之间的相位差是多少,只要对C做FFT变换后,除去C中的50hz分量,就能把B除去?(假设A是不含50hz分量的)
 
to yangzi:

原帖由 aiyou18 于 2007-10-23 17:07 发表
我因为不知道两个波形之间的相位差的问题。所以我才没有在时域里直接相减。
想问一下,在频域里相减需要考虑相位差吗?
yangzi:频域里复数相减和直接在时域里相减是完全一样的。傅立叶变换具有线性性

这个是yangzi你的回答,我觉得不对。

我用MATlab做了下试验。我用一个100hz的正弦波和一个300hz的正弦波叠加:x=sin(2*pi*f1*t)+sin(2*pi*f2*t);我改变这个公式中,两个波的相位差,比如x=sin(2*pi*f1*t+pi/7)+sin(2*pi*f2*t);

然后再对改变相位差前后的波形做FFT变换。我发现幅频曲线是完全一样的,只有100hz与300hz分量。改变相位差不会出现其他的频率分量。

直接在时域里相减需要知道相位差。可是在频域里相减,我可以直接排出掉一个频率分量,然后把它反变换到时域。你看这样思考对吧?
 
--------------------------------------------
yangzj:
回复 aiyou18 的帖子
我说是指两个信号在时域相减后再做FFT,和分别对两个信号做FFT然后再相减,这两个结果是一样的.
频域里相减和直接去掉是两回事.如果能这样做的话那就是一个简单的滤波的问题了.但现在你的问题是信号和噪声是同频率的,你怎么直接去除.
---------------------------------------------
xray:
 
     引用:
另外,请教一个问题,波形A上叠加50hz的正弦波B生成了C。是不是不论A与B之间的相位差是多少,只要对C做FFT变换后,除去C中的50hz分量,就能把B除去?(假设A是不含50hz分量的)
是的。这里,关键是“A是不含50hz分量的”,相当于A和B在频域上是分离的。
-----------------------------------------------
VibrationMaster:
只要波形确定,那么三次谐波就和原始的方波有确定的关系,具体的关系需要自己推。 特别是你的波形带小尖尖,很难在教材中找到现成的答案
0--------------------------------

文章评论0条评论)

登录后参与讨论
我要评论
0
12
关闭 站长推荐上一条 /2 下一条