matlab在振动信号处理中的应用_SVD在信号处理中的应用
写在前面，心烦，还是写点东西转换一下心情。上次说总结一下SVD，今天就写一写，文章公式较多，先写这么多，后续有时间再写，希望你在某天实际应用中，想到我这一篇文章，所以希望你看到后点个赞，谢谢。
不喜欢看理论，matlab代码我也贴出来，可以直接使用。
%%%%%%%%%%%%%%%%%%%%%%
关于SVD分解的基本理论以及其中的物理含义我们不做过多的解读，这里我们只涉及SVD在吸信号处理中的那些奇淫技巧，如果写的不对或者需要补充的可在评论区提出。
关于SVD基本用法如下：
1：线性解调：
关于利用SVD进行相位提取的可以参考我上一篇文章，这里不做介绍。
2：IQ通道效正：
但是在实测过程中，I，Q两路正交通道信号存在“幅度—相位不平衡”的问题，导致上述情况对实测测量造成影响。具体表述为：
对于单频差拍信号
而实际两通道存在“幅度—相位不平衡”。以I路通道为基准，设Q路通道幅度相对增加了
由此合成的差拍信号波形为：
可以看到，由于I,Q两路通道幅相不平衡，合成的差拍信号不仅包含了上边带频率分量，还产生了与上边带频率分量镜像对称的下边带频率分量。
假设得到的回波数据如下所示：
如果,IQ两通道正交，论文中指出，其
存在一种变换T 可以校正 I、Q 通道的相位误差与增益误差，使得
存在幅度相位误差的 I、Q通道的自相关矩阵R ，对采样数据矩阵H 进行奇异值分解得
假设
对于效正后的相关矩阵
可知：
matlab程序如下：
H=[I1;Q1]';
[U,S,V] = svd(H);
T=V*pinv(S(1:2,1:2))*V';
H0=H*T;
x=(H0(:,1)+j*H0(:,2)).';仿真结果如下：（采用实测单手摆动数据）
从结果可以看出，镜像分量有效的得到了抑制。
3：去噪
如果一个平稳信号含有噪声，采用何种方法进行降噪，事实证明SVD也是可以的。
如何将一个一维信号
__
汉克尔矩阵 (Hankel Matrix) 是指每一条副对角线上的元素都相等的方阵，其频繁用于噪声的去除。一维信号构造汉克尔矩阵如下所示：
一般可取，
噪声分离的前提是信号与噪声不相关。如果信号具有高度 self-correlated，其信号的能量就只集中在几个奇异值上，从而达到信号去噪的结果。并且对于复杂的调频信号其奇异值不具有这种特性，故不能去噪。
信号恢复可以从第一列和最后一行得到。
具体的matlab代码如下：
p=round(N/2+1);
q=N+1-p;
H=zeros(p,q);
for ii=1:q
H(:,ii)=s(ii:ii+p-1).';
end
[U,S,V] = svd(H);
S2=zeros(size(S));
for ii=1:10
S2(ii,ii)=S(ii,ii);
end
sr=U*S2*V';
sr2=[sr(1:p,1).' sr(p,2:end)];仿真如下：
原始信号
奇异值
降噪结果：
4：PCA
关于PCA的原理以及实现资料众多，这里就不多赘述，感兴趣的可以自己看看matlab中关于PCA的介绍和使用。
5：多谐波信号频率估计
假设在一个很短的时间内，信号是由一系列任意幅度、频率、衰减、初相的指数函数线性构成的。表达形式如
求解上述模态的频率及相位等信息，可知可以通过FFT进行计算，但是后续的比较可知，有比FFT更为精确的表达。
实际推导过程很复杂，感兴趣的可以看看论文3和论文4，我们这里只是将主要结果，以及如何实现给出。
两篇论文中的方法有点区别，但是很类似，我不知方法是不是等效的，欢迎有关人员继续研究。
论文3中的方法：
主要流程如下：
构造 Hankel矩阵：
row-enhanced data matrix
enhanced Hankel matrix
由于我们处理一维数据，可知
其中：
由于：
，假设模型阶数为p,则
；
在进行特征值分解：
，取其中的特征值为：
，估计频率如下：
代码如下：
%% 第一种
N=400;
p=2;
y=x1;
K=round(N/2+1);
L=N+1-K;
x_n=zeros(L,K);
for ii=1:L
x_n(ii,:)=y(ii+(0:K-1));
end
[s,v,d] = svd(x_n); % X = U*d*V'.
v_sn=v(1:p,1:p);
15. owmeg=s(:,1:p)*sqrt(v_sn);
tao=sqrt(v_sn)*d(:,1:p)';
owega_r1=owmeg(2:N-K+1,:);
owega_rL=owmeg(1:N-K,:);
Ar=inv(owega_rL'*owega_rL)*owega_rL'*owega_r1;
[V,D] = eig(Ar);
for ii=1:p
a(ii)=log(D(ii,ii))*Fs
fi(ii)=angle(D(ii,ii))*Fs/2/pi
end论文4中的方法如下：
第一步，依旧是汉克尔矩阵：
在进行特征值分解：
，取其中的特征值为：
，估计频率如下：
代码如下：
%% 第二种
N=400;
p=2;
y=x1;
K=round(N/2+1);
L=N+1-K-1;
x_n=zeros(L,K);
x_n_2=zeros(L,K);
for ii=1:L
x_n(ii,:)=y(ii+(0:K-1));
x_n_2(ii,:)=y(ii+(0:K-1)+1);
end
[s,v,d] = svd(x_n); % X = U*d*V'.
16. v_sn=v(1:p,1:p);
s1=s(:,1:p);
d1=d(:,1:p);
Ar=v_sn^(-1/2)*s1'*x_n_2*d1*v_sn^(-1/2);
[V,D] = eig(Ar);
for ii=1:p
a(ii)=log(abs(D(ii,ii)))*Fs
fi(ii)=angle(D(ii,ii))*Fs/2/pi
end测试结果：
实验数据：
x1=0.1*exp(j*2*pi*50*t)+0.2*exp(j*2*pi*100*t);结果如下：
方法一：
a =
1.0e+02 *
0.0000 + 6.2832i -0.0000 + 3.1416i
fi =
100.0000 50.0000方法二：
a =
1.0e+02 *
-0.0000 + 6.2832i -0.0000 + 3.1416i
fi =
100.0000 50.0000从结果来看，频率估计是正确的，振幅估计也符合两者的大小关系，所以结果应该是正确的，如果不对，欢迎评论区指出。
参考文献：
[1]Chen S, Yang Y, Wei K, et al. Time-Varying Frequency-Modulated Component Extraction Based on Parameterized Demodulation and Singular Value Decomposition[J]. IEEE Transactions on Instrumentation & Measurement, 2016, 65(2):276-285.
[2]F. Quaiyum, N. Tran, J. E. Piou, O. Kilic and A. E. Fathy, "Noncontact Human Gait Analysis and Limb Joint Tracking Using Doppler Radar," in IEEE Journal of Electromagnetics, RF and Microwaves in Medicine and Biology , vol. 3, no. 1, pp. 61-70, March 2019.
[3] J. E. Piou, K. Naishadham and A. E. Fathy, "A 1-D block processing for non-invasive detection of 2-D cardiac and respiratory rates," 2018 IEEE Radar Conference (RadarConf18) , Oklahoma City, OK, 2018, pp. 0029-0032.
[4] James Hu S L , Yang W L , Li H J . Signal decomposition and reconstruction using complex exponential models[J]. Mechanical Systems and Signal Processing, 2013, 40(2):421-438.