不用傅里叶变换，提取某一频率的幅值和相位

摘要：本文从实际工程问题入手，探寻解决办法，为引入信号正交分解，和广义傅里叶级数做铺垫。

转子做周期性旋转时，不平衡质量所产生的周期性惯性离心力会引起转子产生简谐运动。在工程现场测得某一不平衡转子产生的振动加速度信号(记为 s i g sig sig)如下图所示：

从图中我们可以观测出不平衡分量的频率大约为7.4 Hz, 现在抛出问题，如何提取该信号不平衡成分的幅值和相位？
大家首先想到的肯定是：傅里叶变换

傅里叶变换肯定是可以的，但今天讨论的是，假如没有傅里叶变换，该如何解决？本文会给出求解过程并引出信号处理中的重要理论——信号的正交分解和广义傅里叶级数。求解过程如下：

假设由该不平衡量产生的振动加速度波形(记为 f i t fit fit)可以表示为：
f i t = A ⋅ c o s ( 2 π f t + φ ) fit=A \cdot cos(2\pi ft+\varphi) fit=A⋅cos(2πft+φ)
则该问题转换为在 f = 7.4 f=7.4 f=7.4 Hz时，寻找一对 ( A , φ ) \begin{pmatrix}A ,\varphi\end{pmatrix} (A,φ)使得， f i t fit fit与上述图中的信号 s i g sig sig距离 ϕ \phi ϕ最小：

ϕ ( A , φ ) = ∑ j = 1 N ( f i t [ j ] − s i g [ j ] ) 2 \phi (A, \varphi) = \sum_{j=1}^N(fit[j]-sig[j])^2 ϕ(A,φ)=j=1∑N(fit[j]−sig[j])2
其中 N N N为振动信号 s i g sig sig和 f i t fit fit的点数，为方便计算，对fit进行展开，有
f i t = A ⋅ c o s ( φ ) ⋅ c o s ( 2 π f t ) − A ⋅ s i n ( φ ) ⋅ s i n ( 2 π f t ) fit=A \cdot cos(\varphi) \cdot cos(2\pi ft) - A \cdot sin(\varphi) \cdot sin(2\pi ft) fit=A⋅cos(φ)⋅cos(2πft)−A⋅sin(φ)⋅sin(2πft)

令 a = A ⋅ c o s ( φ ) a = A \cdot cos(\varphi) a=A⋅cos(φ) , b = − A ⋅ s i n ( φ ) b = - A \cdot sin(\varphi) b=−A⋅sin(φ), x = c o s ( 2 π f t ) x = cos(2\pi ft) x=cos(2πft), y = s i n ( 2 π f t ) y= sin(2\pi ft) y=sin(2πft), 则有

f i t = a ⋅ x + b ⋅ y fit= a \cdot x + b \cdot y fit=a⋅x+b⋅y
因为 f f f 是确定的，采样时刻序列 t t t 也是确定的，所以 x x x 和 y y y 都是常量序列。则此时 f i t fit fit 与 s i g sig sig 的距离 $ \phi $ 可以表示为以 a a a， b b b 为变量的表达式：

ϕ ( a , b ) = ∑ j = 1 N ( f i t [ j ] − s i g [ j ] ) 2 \phi (a, b) = \sum_{j=1}^N(fit[j]-sig[j])^2 ϕ(a,b)=j=1∑N(fit[j]−sig[j])2

将 f i t = a ⋅ x + b ⋅ y fit= a \cdot x + b \cdot y fit=a⋅x+b⋅y 代入上式，有：
ϕ ( a , b ) = ∑ j = 1 N ( a ⋅ x [ j ] + b ⋅ y [ j ] − s i g [ j ] ) 2 \phi (a, b) = \sum_{j=1}^N(a \cdot x[j] + b \cdot y[j] -sig[j])^2 ϕ(a,b)=j=1∑N(a⋅x[j]+b⋅y[j]−sig[j])2
将上式完全展开，有：

ϕ ( a , b ) = ∑ j = 1 N ( a 2 ⋅ x 2 [ j ] + b 2 ⋅ y 2 [ j ] + s i g 2 [ j ] + 2 a b x [ j ] y [ j ] − 2 a ⋅ x [ j ] s i g [ j ] − 2 b ⋅ y [ j ] s i g [ j ] ) \phi (a, b) = \sum_{j=1}^N(a^2 \cdot x^2[j] + b^2 \cdot y^2[j] + sig^2[j] + 2abx[j]y[j] - 2a \cdot x[j]sig[j] - 2b \cdot y[j] sig[j]) ϕ(a,b)=j=1∑N(a2⋅x2[j]+b2⋅y2[j]+sig2[j]+2abx[j]y[j]−2a⋅x[j]sig[j]−2b⋅y[j]sig[j])

根据连续函数极值定理，在 ϕ \phi ϕ 最小时，有：
∂ ϕ ∂ a = ∂ ϕ ∂ b = 0 \frac{ \partial \phi }{ \partial a } = \frac{ \partial \phi }{ \partial b } = 0 ∂a∂ϕ=∂b∂ϕ=0

整理，得到方程组：

{ ∑ j = 1 N x 2 [ j ] ⋅ a + ∑ j = 1 N x [ j ] y [ j ] ⋅ b = ∑ j = 1 N s i g [ j ] x [ j ] ∑ j = 1 N y [ j ] x [ j ] ⋅ a + ∑ j = 1 N y 2 [ j ] ⋅ b = ∑ j = 1 N s i g [ j ] y [ j ] \begin{cases} \sum_{j=1}^Nx^2[j]\cdot a +\sum_{j=1}^Nx[j]y[j] \cdot b=\sum_{j=1}^Nsig[j]x[j] \\ \sum_{j=1}^Ny[j]x[j] \cdot a +\sum_{j=1}^Ny^2[j]\cdot b=\sum_{j=1}^Nsig[j]y[j] \end{cases} {∑j=1Nx2[j]⋅a+∑j=1Nx[j]y[j]⋅b=∑j=1Nsig[j]x[j]∑j=1Ny[j]x[j]⋅a+∑j=1Ny2[j]⋅b=∑j=1Nsig[j]y[j]

令：

X = [ x [ 1 ] y [ 1 ] x [ 2 ] y [ 2 ] . . . . . . x [ N ] y [ N ] ] , Y = [ s i g [ 1 ] s i g [ 2 ] . . . s i g [ N ] ] , α = [ a , b ] X = \begin{bmatrix}x[1] & y[1]\\x[2] & y[2]\\... & ...\\x[N] & y[N] \end{bmatrix},Y = \begin{bmatrix} sig[1]\\sig[2]\\...\\sig[N]\end{bmatrix}, \alpha = [a, b] X=⎣⎢⎢⎡x[1]x[2]...x[N]y[1]y[2]...y[N]⎦⎥⎥⎤,Y=⎣⎢⎢⎡sig[1]sig[2]...sig[N]⎦⎥⎥⎤,α=[a,b]

则上述方程组可表示为：

X T X α = X T Y X^TX \alpha = X^TY XTXα=XTY

进而求得：
α = ( X T X ) − 1 X T Y \alpha = (X^TX)^{-1}X^TY α=(XTX)−1XTY

求得 α \alpha α 后，便得到了 a , b a, b a,b 的值，进而通过下面的公式得到 A , φ A , \varphi A,φ 的值：
A = a 2 + b 2 , φ = − a r c t a n b a A = \sqrt{a^2 + b^2}, \varphi = -arctan\frac{b}{a} A=a2+b2 ,φ=−arctanab
至此，问题得到解决，上述方法成功提取到了该信号不平衡成分的幅值和相位。效果见下图所示，核心代码在文末给出。

这时，有人会问了，你是怎么观测到不平衡分量的频率大约为7.4 Hz的？我看着怎么有点像7.3Hz呢？其实相差不大，下图是本文分析过程分别在f为7.4, 7.3, 7.2Hz时的结果对比图。

那么如何理解这种结果呢？

其实，当f等于7.4Hz时，所提取的信号成分是原始信号在"7.4Hz方向"上的投影。7.4Hz方向和7.3Hz方向以及7.2Hz方向的"夹角"很小，所以投影的结果差距也不是很大，因此在精度要求不是特别苛刻时，可以使用肉眼观测的频率(7.4, 7.3, 7.2都可以)进行提取。

通过上面的论述可以知道，信号可以在任意多个方向上进行投影，在很多时候，我们需要的可能是信号在多个方向上的投影，那么有没有一种通用的大家都认可的投影方式呢？

答案是肯定的，参考二维平面中向量的表示方法或许可以得到一些启示：二维平面中的向量可以用一对坐标来表示，例如（ 3 , 4 ）（3,4）（3,4）代表该向量在x轴方向上投影的长度为3，在y轴方向上投影的长度为4. 那么该向量可以分解为x轴方向和y轴方向的两个向量（ 3 , 0 ）（3,0）（3,0）和（ 0 , 4 ）（0,4）（0,4），在中学学习受力分析时我们也是这么干的。
那么为什么沿着x轴方向和y轴方向呢？因为这两个方向是正交的。而正交分解得到的分量具有很多美好的性质，在下一篇文章中，会深入浅出地介绍信号的正交分解和广义傅里叶级数。
微信公众号“振动信号研究所”会在第一时间发布新作，关注后回复“unbalance”获取本文所用波形数据。

import numpy as np
import matplotlib.pyplot as pltsig = np.loadtxt(r".\unbalanced.txt")
Fs = 8192
t = np.arange(len(sig))/Fs
f = np.array([7.4 for _ in t])  # 欲提取的频率成分
x = np.cos(f*t*2*np.pi)
y = np.sin(f*t*2*np.pi)X = np.mat(np.asarray([x, y]).T)
Y = np.mat(sig).T
xtx = X.T * Xalpha = xtx.I * (X.T * Y)a = alpha.tolist()[0][0]
b = alpha.tolist()[1][0]
A = np.sqrt(a**2 + b**2)
phi = - np.arctan(b/a) if - np.arctan(b/a) < 0 else -np.pi - np.arctan(b/a)
fit = A*np.cos(2*np.pi*f*t + phi)
plt.figure()
plt.plot(t, sig, label="sig")
plt.plot(t, fit, label="fit")
plt.title(f"fit: {np.round(A, 3)}*np.cos(2*np.pi*{f[0]}*t+({np.round(phi, 3)}))")
plt.legend(loc="upper right")
plt.show()

不用傅里叶变换，提取某一频率的幅值和相位相关推荐

下垂控制的两电平三相桥式逆变器，可根据负载的变化改变三相输出电压的频率和幅值
下垂控制的两电平三相桥式逆变器,可根据负载的变化改变三相输出电压的频率和幅值. ABC到dq解耦后,采用电压电流双闭环控制,下垂控制给定坐标系变换的频率和幅值. 几年前做的仿真,支持MATLAB201 ...
使用数字示波器DS6104测量交流信号的幅值和相位
01简介使用普通的万用表测量交流信号的时候,通常会遇到万用表的频率响应的问题.使用可以联网的示波器可以获得它采集到的数据,进而可以计算出所测量的交流信号的有效值和相位. 这里通过实验来确定使用示 ...
matlab复数的相位,复数的幅值和相位
[i]); x[i]=x[i]*180/PI-90; printf("第%d次谐波的相位为为%f \\n", i , x[i]) ; } printf("第%d次谐波的幅 ...
matlab 求电流幅值,输出信号的幅值与相位.ppt
输出信号的幅值与相位 1.什么叫频率分析?2.频率分析的问题引入?3.频率分析的仿真.4.频率分析的计算方法.5.频率分析的作用与意义6.幅频特性与相频特性王选择 1 频率分析的定义 1 频率分析的 ...
matlab频率和幅值图,请教！傅里叶变换频率和幅值对不上！
拜托各位大哥帮忙看看,这个程序运行了之后幅值和应有频率对不上啊?自己找不到问题.而且相位图也不对. clc; close all; clear; A1=2.0; %频率F1信号的扰动幅度 A2=1.9 ...
基于STM32F407实现快速傅里叶变化（FFT），计算指定频率的幅值
本人的课题是关于EIT采集系统,简单的说就是往人体注入特定频率的恒流源,再采集电压信号,通过分析电阻抗分布进行成像.采集的电压信号是需要进行FFT处理,只保留注入频率的信号成分.本文主要介绍如何在ST ...
傅里叶变换的理解----计算幅值和相位
先推荐表格文章如果看了此文你还不懂傅里叶变换,那就过来掐死我吧[完整版] 傅里叶原理表明:任何连续测量的时序或信号,都可以表示为不同频率的正弦波或余弦信号的无限叠加. 1.FT的理论就会告诉你可以 ...
FFT之频率与幅值的确定
FFT也就是快速傅里叶变换.经过快速傅里叶变换后会得到一串复数. 下面要讲两个问题:1.如何获取频率:2.如何获取幅值傅里叶变换并没对频率进行任何计算,频率只与采样率和进行傅里叶变换的点数相关. F ...
怎么检测声音频率和幅值_【电缆小课堂】国网天津电缆公司电缆带电检测团队介绍及典型案例分析...
●团队介绍● 国网天津市电力公司电缆分公司电缆运检北部中心电缆试验班成立于2019年9月,主要从事国网天津电缆公司辖区内高压电缆带电检测工作.班组现有员工11名,其中硕士8人,占总人数72.8%,35 ...

不用傅里叶变换，提取某一频率的幅值和相位

不用傅里叶变换，提取某一频率的幅值和相位

不用傅里叶变换，提取某一频率的幅值和相位相关推荐

最新文章

热门文章