精华内容
下载资源
问答
  • 行业分类-电信-同步分离信号与噪声的信号分析装置及其分析方法.rar
  • 内容简介: 脑电信号分析已经在脑科学研究中占据了越来越重要的地位...《脑电信号分析方法及其应用》可供生物医学工程中脑信号处理方面的研究人员、大中专院校的相关专业的研究生,以及医院脑电图室的医务工作者参考。
  • 本文将讨论开关稳压器PSRR(电源抑制比,其对输入噪声抑制很重要)以及信号分析方法。开关纹波噪声本部分依据基波和谐波理论介绍降压转换器输出纹波计算公式。根据开关稳压器拓扑结构和基本操作,纹波始终是开关稳压...
  • 内燃机的非平稳信号分析方法及其噪声源小波识别技术的研究
  • 浅海中的宽带水声信号传播呈现出频散的特点,通过高分辨率的时频分析方法可以刻画频散曲线。通过数值仿真和实验数据处理,对比分析几类常用的时频分析方法在提取宽带声信号频散曲线方面的性能。结果表明:STFT时频...
  • 随机信号实验--微弱信号的检测提取及分析方法 基于多重自相关
  • 阐述了心音信号的产生机理、组成成分以及人体微弱心音信号检测的关键技术,介绍了心音信号处理技术在心血管疾病无创诊断中的意义,结合应用分析了经典心音信号分析方法的局限性,对现代心音分析中的常用的时频分析...
  • 行业分类-电子电器-无线电信号分析方法及装置.zip
  • MATLAB不同时频信号处理方法介绍效果对比

    千次阅读 多人点赞 2018-10-25 17:29:40
    目录MATLAB不同时频信号处理方法介绍效果对比信号的频域分析Fourier(FFT)变换基本原理MATLAB 中信号FFT的基本处理方法MATLAB 中不同参数设置对FFT结果的影响信号的功率谱分析的基本原理MATLAB中信号功率谱分析的...

    本文欢迎非商业目的的学习分享转载,转载请附上原文链接及作者ID

    本文为作者自身的一个学习总结,大部分内容在相关教材上也可以找到,有空的也会不定期更新。本身也在学习的过程中,出现错误在所难免,欢迎大家在留言区指正。看到的话我都会及时回复的。
    文章中使用的数据有需要的话可以从此处进行下载

    MATLAB不同时频信号处理方法介绍及效果对比

    从大二到现在接触傅里叶变换已经好几年了,对傅里叶变换的基本概念与方法和使用的注意事项也都有所了解,但动手方面一直较为欠缺,最近刚好又重新在做相关方面的可以,顺手通过Matlab把整个信号的FFT和相关的时频概念都串起来复习一下,方便以后查找。

    不同于一般教程中使用生成的函数进行分析,本文直接使用一段实际采集的电机径向的振动信号进行分析处理,更真实地符合现实情况

    1. 信号的频域分析

    1.1 Fourier(FFT)变换基本原理

    傅里叶变换的基本公式基本上可以在任何一本有关信号处理的大学教材或者相关的培训PPT上找到,在这里不再赘述。而关于傅里叶变换究竟是什么,网上有很多通俗化的解释,看到的比较好的文章有知乎的一篇专栏:如果看了这篇文章你还不懂傅里叶变换,那就过来掐死我吧,另外b站上也有一个国外大神3Blue1Brown制作的一个可视化效果非常出众的短片:形象展示傅里叶变换
    傅里叶变换的可视化表示:Alt

    简单来说就是任何一个时域信号都可以用一系列正弦函数sin(nx)的和来表示,并且使用的正弦函数项越多,所得到的“模拟信号”也就更接近真实信号。信号的傅里叶变换实际上就是求这一系列正弦函数系数的过程。

    至于为什么可以对一个信号进行傅里叶变换和为什么能这么对信号进行变换这个问题我一直到研一上完“矩阵分析”,和“数字信号处理”这么课之后才有了比较清楚的认识。无论是傅里叶变换,短时傅里叶变换,还是后面会谈到的小波变换,经验模态分解。其本质都是时域信号在线性空间正交基上的分解与重构,而信号分析的重点则在于如何针对不同的信号来选择相应合适的线性空间正交基,或者说各类变换的“核函数”

    至于为什么学的时候教的都是傅里叶变换,而实际用的函数都是FFT(快速傅里叶变换),这是

    MATLAB 中信号FFT的基本处理方法

    在MATLAB中进行FFT的函数为fft()

    y=fft(x,n)
    

    其中x为时域信号,参数n表示执行n点fft变换
    下面举一个最简单的例子

    clear;clc;
    tic
    load('vibdata.mat')
    fs=5120;%实际采样频率
    nfft=2^12;%fft变换点数4096 
    x=fft_data(1:2^12);%取1秒不到的时域数据
    fy=abs(fft(x,nfft));
    fx=(1:nfft/2)*fs/nfft;
    t=0:1/fs:nfft/fs-1/fs;
    figure()
    subplot(121)
    plot(t,x);
    xlabel('t/s');ylabel('幅值');title('原始时域信号')axis([0 nfft/fs -1 1])
    subplot(122)
    plot(fx,fy(1:nfft/2)*2/nfft);
    xlabel('频率/Hz');ylabel('幅值');title('FFT频域信号');axis([0 fs/2 0 0.1])
    toc
    %时间已过 0.098931 秒
    

    MATLAB 中不同参数设置对FFT结果的影响

    1. 采样点数的影响

    根据香农采样定律,若已知信号的最高频率为fc,为防止混叠,则选定抽样频率fs应该满足fs≥2fc,同时根据选定的频率分辨率f0,则确定进行FFT变换所需要的点数N=fs/f0=Tfs

    在采集信号的过程中,为了防止混叠,首先应选定抽样频率fc,但同时该频率不建议选得太高,否则为了实现理想的频域分辨率,N会变得很大。极大地增加FFT过程中的计算量。

    现在对上述FFT程序做适当更改,比较当数据点数分别等于2048,4096和8192时FFT的变换结果。

    clear;clc;
    tic
    load('vibdata.mat')
    fs=5120;%实际采样频率
    nfft=[2^11 2^12 2^13];%fft变换点数2048,4096,8192
    ndata=[2^11 2^12 2^13];%数据量,取1秒,2秒,3秒左右的数据
    x1=fft_data(2^12:ndata(1)+4095);fy1=abs(fft(x1,nfft(1)));fx1=[(1:nfft(1)/2)*fs/nfft(1)];
    x2=fft_data(2^12:ndata(2)+4095);fy2=abs(fft(x2,nfft(2)));fx2=[(1:nfft(2)/2)*fs/nfft(2)];
    x3=fft_data(2^12:ndata(3)+4095);fy3=abs(fft(x3,nfft(3)));fx3=[(1:nfft(3)/2)*fs/nfft(3)];
    figure()
    plot(fx1,fy1(1:nfft(1)/2)*2/nfft(1),'-r',...
        fx2,fy2(1:nfft(2)/2)*2/nfft(2),'-b',...
        fx3,fy3(1:nfft(3)/2)*2/nfft(3),'-g');
    xlabel('频率/Hz');ylabel('幅值');title('FFT频域信号');axis([0 fs/2 0 0.1]) 
    legend('2048','4096','8192');
    toc
    


    从结果中切实可以比较清楚的看出频率分辨率随着处理点数的增加而增加了,原始信号为30Hz的振动信号,当采样频率增加之后,原本被淹没掉的谐波分量也显示了出来。
    在这里插入图片描述
    但是需要注意的是,当采样时间较长时,系统状态可能发生改变,虽然频谱“分辨率”提高了,但是又重新引入了不少杂质分量,因此不一定能够提高频率的分辨率

    2. fft变换点数nfft对频率分辨率的影响

    在MATLAB中fft函数的基本形式是y=fft(x,n),其中n为fft变换的点数,一般来说为了提高fft计算的效率,n都会选择一个2的整数次幂的数值,这里假定该数值为nfft,信号x时域数据中点点数为N。一般来说为了不浪费时域在MATLAB中fft函数的基本形式是y=fft(x,n),其中n为fft变换的点数,一般来说为了提高fft计算的效率,n都会选择一个2的整数次幂的数值,这里假定该数值为nfft,信号x时域数据中点点数为N。一般来说为了不浪费时域中的信息,nfft≥N,多出来的点数用0补足,这个操作我们称之为补零(zero padding)

    在前面有说过频谱的分辨率f0=fs/nfft,nfft为进行fft变换的点数,当我们进行补零操作后,参与fft变换的点数有N变为nfft,变大了,看起来频谱的分辨率f0增加了,我们在没有增加数据量的情况下提高了频谱的分辨率,上述说法是不准确的

    这里有两个概念需要说明,一个是物理分辨率,一个是计算分辨率

    物理分辨率:f0=1/T Hz,只和采样时间有关
    计算分辨率:f0=fs/nfft Hz,和fft变换点数有关

    这里借用一篇文章中的几幅图片对这个概念进行一下说明。下面第一幅图是原始数据图像,第二幅图是补零后的数据图像。
    在这里插入图片描述在这里插入图片描述 对比两幅图片,其实可以发现补零后的fft的分析对象发生了改变。我们知道实际信号一般是无限长序列,而FFT针对的序列总是有限长度的,因此,DFT假设实际信号为对输入的有限长序列信号的周期延拓。将上面两幅图分别进行周期延拓之后,可以得到下面的周期序列。
    在这里插入图片描述
    在这里插入图片描述
    对比上面两幅图,其实可以看到对于同样长度的分析数据,通过补零得到的数据可以看做对原始数据加窗后得到的数据。而对数据加窗就会在频域上产生频谱泄露,这也就解释了为什么补零能提高频域的分辨率,但却无法提供更多的信息

    下面在举一个例子来说明这个问题,以256Hz的采样率对70Hz的正弦波采集1s的数据,并对该段数据做FFT变换,可以得到下列图像,可以看到FFT准确地描述了原始信号的频谱分量。
    FFT变换图像
    当我们队原始信号进行补零操作后,观察得到的FFT图像如下。可以非常明显的看到频谱泄露的现象。当我们把70Hz附近的图像放大之后,可以看到原本的谱线变成了类似于sinc(f)的形式。这也就印证了前面印证前面加窗的说法。因为时域上的加窗等于频域上的卷积,而矩形窗FFT的结果就是一个sinc()函数。


    对于上述问题有一种解决方案是用本身的数据代替零点来补足缺少的问题,这种方法也存在一些问题,这里不做讨论,而是介绍另外一种加窗的操作方法。通过对补零后的数据添加类似hann(hanning窗)等窗函数,来减小频率泄露的问题。这样做的原因是时域加窗等于频域上的卷积,而合适的窗函数在频域上的能量较为集中,频谱泄露的情况也能有所减轻。但加窗又会出现一些其他的问题,可以参考我的另外一边文章
    FFT快速傅里叶变换中的误差来源》


    最后我们来看一下不同补零形式所体现的效果。从上到下分别是未补零,直接补零和加窗补零。

    最后还是通过MATLAB来实践一下不同nfft点数的FFT处理效果。

    clear;clc;
    tic
    load('vibdata.mat')
    fs=5120;%实际采样频率
    nfft=[2^11 2^12];%fft变换点数2048,4096,8192
    ndata=[2^11 2^12];%数据量,取1秒,2秒,3秒左右的数据,4秒左右的数据
    x1=fft_data(1:ndata(1));fy1=abs(fft(x1,nfft(1)));fx1=[(1:nfft(1)/2)*fs/nfft(1)];
    x2=fft_data(1:ndata(1));fy2=abs(fft(x2,nfft(2)));fx2=[(1:nfft(2)/2)*fs/nfft(2)];
    x3=fft_data(1:ndata(1)).*hann(2^11);fy3=abs(fft(x3,nfft(2)));fx3=[(1:nfft(2)/2)*fs/nfft(2)];
    x4=fft_data(1:ndata(2));fy4=abs(fft(x4,nfft(2)));fx4=[(1:nfft(2)/2)*fs/nfft(2)];
    figure(1)
    plot(fx1,fy1(1:nfft(1)/2)*2/nfft(1),'-r',...
        fx2,fy2(1:nfft(2)/2)*2/nfft(2),'-b',...
        fx3,fy3(1:nfft(2)/2)*2/nfft(2),'-g',...
        fx4,fy4(1:nfft(2)/2)*2/nfft(2),'-c');
    xlabel('频率/Hz');ylabel('幅值');title('FFT频域信号');axis([0 fs/2 0 0.1]) 
    legend('未补零','直接补零','加窗补零','增加采样点');
    figure(2)
    t2=0:1/fs:4095/fs;
    plot(t2,x4,':b',t2,[x3;zeros(2048,1)],'-r',t2,[hann(2048);zeros(2048,1)],'-g');
    xlabel('时间/s');ylabel('幅值');title('时域信号');
    %axis([0 fs/2 0 0.1]) 
    legend('原始两倍信号','加窗后的补零信号','hanning窗');
    toc
    

    在这里插入图片描述
    在这里插入图片描述
    通过观察上述fft处理后的图像可以发现,与原始信号相比,直接补零虽然看似提高了频域的分辨率,但也造成了较明显的频谱泄露现象。而如果以增加FFT点数后的FFT变换效果为参考,加窗补零这种方法不仅在一定程度上提高了频谱的分辨率,而且与直接补零相比还叫好地抑制了频谱泄露的现象。

    说了这么久补零的局限性,最后简单说一下补零的优点:除了上述提到的(1)使得数据长为2的整数次幂以提高FFT的算法外(2)克服栅栏效应,是谱的外观得到平滑,提高谱的分析进度。

    这里盗一张我们老师的课程PPT来进行说明,虽然补零并不能真正地提高频率分辨率,但是在非正周期采样的过程中,但是通过合理的补零可以对连续谱的谱峰值进行有效地逼近。
    在这里插入图片描述

    3. 窗函数对FFT变换效果的影响

    在这里插入图片描述
    这一小节主要对窗函数对作用进行分析,在进行数字信号处理时,不可能对无限长(完整)信号进行测量和运算,只能取有限的时间片段进行分析,这个过程称为信号截断,时域上的截断等于给原始信号乘上一个窗函数。而时域上两个函数相乘,在频域就是其频谱的卷积。由于窗函数不可能无限宽,其频谱不可能为冲激励函数,信号的频谱与窗函数的卷积必然产生拖尾现象,造成频谱泄露。原始谱线的能量因为与窗函数的卷积被分散到了两个较宽的频带当中,造成了能量泄露

    频谱泄漏是由FFT算法中的一个假设导致的,即持续精确地重复时间记录,且时间记录中包含的信号在对应时间记录长度的间隔内呈周期性。若时间记录的周期数为非整数, 便违反该假设,导致频谱泄漏。另一种看法是,信号的非整数周期频率分量未与频谱频率线之一精确吻合

    您仅可在两种情况下保证获得整数周期。第一种情况是,您对测量的信号进行同步采样,因此可按需获取整数周期。

    另一种情况是,您捕获的瞬时信号可完全融入时间记录。但是,多数情况下,您测量的未知信号是平稳的,即该信号在采集前、中、后都存在。这时,您就无法保证采样的是整数周期。频谱泄漏对测量造成干扰,来自给定频率分量的能量将分散至相邻的频率线或仓。您可使用窗口,将在非整数周期内进行FFT产生的效果最小化。

    以矩形窗为例,有限宽度的矩形窗在频谱中的图像为sinc()函数。主瓣的宽度会随着N采样点数的增加而减小,从而降低频域卷积过程中造成的能量泄露。

    N = 8; xn = ones(1,N);
    Nfft = 8192 %较密的谱线近似连续谱
    w = [-pi:2*pi/Nfft:(pi-2*pi/Nfft)];
    Xk = fft(xn,Nfft);
    figure
    plot(w,abs(fftshift(Xk)))
    

    在这里插入图片描述
    因为在进行FFT的过程中,时域的截断是必然的,因此频域泄露也是必然的。因此只能想办法来减少截断过程中造成的频域泄露现象。主要有以下两种方法:(1)加大窗口宽度,既增加采样点数N;(2)采用适当的窗函数进行截断;

    那什么是适当的窗函数呢?,一般认为窗函数的主瓣宽度应该尽可能的小,以提高频率分辨率,同时主瓣和第一旁瓣的幅值差值应该尽可能的大,并且旁瓣的衰减越快越好,从而减小频谱泄露提高频谱幅值精度

    但“鱼和熊掌不可兼得”的朴素真理在这里同样适用,主瓣宽度与旁瓣的幅度的衰减速度相互制约。主瓣宽度越小,频率分辨率越高,但旁瓣幅值较大,频谱泄露增加导致频谱精度降低,反之亦然。

    如果仅仅需要求精确获得信号频率,而不考虑幅值精度,例如在进行固有频率测量的时候,就可以选择主瓣较窄的矩形窗。而如果在分析带有较强的干扰噪声的窄带信号的过程中,则应该选择旁瓣幅度较小的汉宁窗,三角窗等。在实际选择窗函数的过程中,需要对各个因素进行权衡处理。下图为各个窗函数的时域和频域波形,主瓣越大的窗,旁瓣衰减越快。
    在这里插入图片描述

    clear;clc;
    tic
    load('vibdata.mat')
    fs=5120;%实际采样频率
    nfft=[2^12 2^13];%fft变换点数4096,8192
    ndata=[2^12 2^13];%数据量,取1秒,2秒,3秒左右的数据,4秒左右的数据
    x1=fft_data(1:ndata(1));fy1=abs(fft(x1,nfft(1)));fx1=[(1:nfft(1)/2)*fs/nfft(1)];
    x2=fft_data(1:ndata(1)).*hann(ndata(1));fy2=abs(fft(x2,nfft(1)));fx2=[(1:nfft(1)/2)*fs/nfft(1)];
    x3=fft_data(1:ndata(1)).*blackman(ndata(1));fy3=abs(fft(x3,nfft(1)));fx3=[(1:nfft(1)/2)*fs/nfft(1)];
    x4=fft_data(1:ndata(1)).*gausswin(ndata(1));fy4=abs(fft(x4,nfft(1)));fx4=[(1:nfft(1)/2)*fs/nfft(1)];
    x5=fft_data(1:ndata(2));fy5=abs(fft(x5,nfft(2)));fx5=[(1:nfft(2)/2)*fs/nfft(2)];
    figure(1)
    plot(fx1,fy1(1:nfft(1)/2)*2/nfft(1),'-r',...
        fx2,fy2(1:nfft(1)/2)*2/nfft(1),'-b',...
        fx3,fy3(1:nfft(1)/2)*2/nfft(1),'-g',...
        fx4,fy4(1:nfft(1)/2)*2/nfft(1),'-m',...
        fx5,fy5(1:nfft(2)/2)*2/nfft(2),'-c');
    xlabel('频率/Hz');ylabel('幅值');title('FFT频域信号');axis([0 fs/2 0 0.1]) 
    legend('矩形窗','Hanning Window','Blackman Window','Gassian Window','两倍矩形窗');
    toc
        
    

    通过对实际信号进行分析可以发现对于实际信号来说,与不加窗(也可以说是加矩形窗)的频谱相比,加了窗的频谱相对而言较为接近,由于主板较宽,频谱泄露因此频谱较为平缓,部分谐波分量被淹没。
    在这里插入图片描述

    4. 频谱细化技术(Zoom-FFT)

    频谱细化指的是在频谱分析的过程中增加频谱中某个频率分辨率的方法,如果我们希望将频段[f1 ~ f2]内的频率分辨率提高K倍,最简单的方法就是在采样频率不变的情况下,将采样点数从N点增加到KN点,但这么做的过程中也同时将许多不感兴趣的频率分辨率也提高了,造成了算力的浪费。
    在这里插入图片描述
    重采样的基本原理如上,在进行细化分析的过程中,必须保证原始信号有足够的长度,如果想把频段的分辨率提高K倍,就得保证原信号的采样长度为KN,这样在每隔K点重采样的过程后,得到新长度为N的复序列频谱细化技术实际上一种降低计算量的方法,而不是直接使用原始数据来提高频率分辨率的方法

    下面通过代码具体分析频谱细化的具体作用:

    此部分暂时无法实现,后续进行内容补充

    信号的功率谱分析的基本原理

    有很长的一段时间内我都混淆了信号的FFT谱与功率谱之间的概念,关于信号的FFT和功率谱可以参考我博客中的两篇文章:《功率谱和频谱的区别、联系》《功率谱估计》。对功率谱与FFT的基本概念有了比较准地说明。

    从个人的角度理解,从名字就能看出这两者的主要区别:“FFT-快速傅里叶变换”,“功率谱估计”。FFT本质上是一个对信号分解的过程,以一种工具。而“功率谱估计”则更侧重估计二字,功率谱估计的出发点就是将分析的对象看做是随机信号,希望通过各种方法来准确地通过采样信号准确地估计出真实信号在频谱上的功率分布,而FFT则是假定分析的对象是一个“确定性的周期信号”来直接对其分解,求解相关的系数。

    而在FFT分析的过程中,为了准确描述信号频谱而采取的类似加窗的操作,个人也有几分估计的味道在里面,并且这种“估计”的操作在功率谱估计的过程中也经常能够看到。同时功率谱估计过程中使用的类似平均等手段个人感觉也可以用在FFT上,但目前这些手段主要还是运用在功率谱估计上,暂时不清楚为何,姑且追寻主流使用这种方法吧

    对于随机信号x(n)来说,实际工作中我们只能得到有限长度N的数据,要由这N个数据来估计信号的均值,方差,自相关函数,功率谱等其他感兴趣的参数。作为一个估计量,我们就需要评价这个估计量的质量,一般的评价手段有偏差,方差,均方差等几个角度。

    经典的非参数功率谱估计方法无法同时在偏差,方差等几个角度达到非常理想的估计效果,这也是现代的功率谱估计方法引入参数法功率估计的原因从这个角度也可以看出功率谱估计和FFT之间的一个显著的区别,我们对信号进行FFT后一般情况下并不会非常在意FFT是否非常准确地还原了原始信号,而是主要关注了FFT频谱中的主要分量,而在使用功率谱估计的过程中,则会较多地关注相关方法的“准确程度”

    关于功率谱估计的质量的具体数学推导过程可以参考比较早的一篇博文《功率谱估计》,在这里不详细描述,只是简单提一下相关的结果。

    MATLAB中信号功率谱分析的基本方法

    经典非参数法功率谱估计

    1. 周期图法

    定义:计算随机信号x(n)的N点观察数据的DFT,然后取其幅值的平方并除以N作为对x(n)真实功率谱的估计

    对于这种方法来说,周期图法估计的方差不随采样点数N的增加而减小。不是功率谱的一致估计,而当N趋向无穷时,变差逐渐趋近于零,是真实值的渐进无偏估计。

    %直接编程
    clear;clc;
    tic
    load('vibdata.mat')
    fs=5120;%实际采样频率
    Nfft=8192;
    x1=fft_data(1:Nfft);
    fy1=abs(fft(x1))/length(x1);
    
    Pxx=2*fy1.^2;
    index=0:round(Nfft/2-1);
    w=index*fs/Nfft;
    figure
    plot(w,10*log(Pxx(index+1)))
     xlabel('频率/Hz');ylabel('功率密度/dB');title('信号的功率谱估计');axis([0 2560 -300 0]); 
    toc
    

    在这里插入图片描述

    %使用MATLAB官方函数
    clear;clc;
    tic
    load('vibdata.mat')
    fs=5120;%实际采样频率
    Nfft=4096;
    x1=fft_data(1:4096);
    [Pww,F,pxxc]=periodogram(x1,[],[],fs,'ConfidenceLevel',0.95);
    figure
    plot(F,10*log(Pww),F,10*log(pxxc),'-.')
     xlabel('频率/Hz');ylabel('功率密度dB/Hz');title('95%置信区间的信号的功率谱估计');axis([0 fs/2 -300 0]); 
    toc
    

    在这里插入图片描述

    这里有一部分问题就是,在计算的过程中,仔细对比两幅图像,可以发现虽然整体形状一致,但手动编程与直接使用MATLAB函数相比,产生PSD整体偏大,该问题暂时自己没法解释,进行搁置

    2. 基于维纳辛钦定律(Wiener-Khintchine theorem)的相关法

    定义:信号的功率谱估计等于该信号自相关函数的离散DTFT(离散时间傅里叶变换) 同时对于平稳随机过程来说,信号的自相关函数就是功率谱密度函数的逆傅里叶变换,而对于非平稳随机过程来说,功率谱密度的逆变换的结果是自想换函数在时域上的平均。
    在这里插入图片描述
    对于相关法来说,当时延m取值的最大值为N-1时,相关法和周期法等效。随着m取值的增大,相关法功率谱估计质量的方差特性和周期图法功率谱估计会具有一样的缺点。

    3. 平滑周期法(Blackman-Tukey法)

    为了降低估计的方差,可以对自相关函数的估计值进行加窗后再进行FFT变换,可以证明这种方法的方差要小于原始的相关法,但这种方法同时也降低了频率的分辨率。

    4. 平均周期法 (Welch-Bartlett法)

    这种方法是将随机序列分为L段,分别求每一段的周期图,最后再进行平均。可以证明当L趋于无穷时,方差趋向于0,是一致估计。但平均图方差减小的代价就是偏差增大。这种方法减小方差的同时也降低了谱的分辨率。主要原因是分段即加窗,窗越短,主瓣宽度越大。

    5. 重叠平均周期法 (Welch法)

    与平均周期法相比,这种方法将各段数据进行了一定的重叠
    在MATLAB中,Bartlett法和Welch法都可以用pwelch函数进行操作。关于pwelch函数的使用方法既可以直接参考MATLAB官方文档,也可以参考我博客中的文章《MATLAB中应用pwelch()函数实现功率谱估计》

    对于几种方法,如果采用一个已知信号进行对比,可以发现:

    1. 周期图法谱估计曲线的波动很大,即估计的方差较大
    2. Welch 法谱估计曲线较为光滑,方差减小,但分辨率降低
    3. 对于Welch 法,当数据分段增加时,各段数据较短时,谱的分辨率明显下降,而谱估计曲线较为光滑,方差较小;反之当数据分段数减小,各段数据较长时,谱分辨率明显提高,而谱估计曲线波动较大,方差较大
      下面采用以一段已知信号来演示各种功率谱估计方法和FFT变换之间的区别:
    clear;clc;
    tic
    fs=1024;
    T=4;
    t=(0:T*fs-1)/fs;
    x=0.1*cos(2*pi*4*t)+0.1*cos(2*pi*4.5*t)+rand(1,size(t,2));
    Nfft=2^nextpow2(size(x,2));
    %FFT变换
    ffy1=abs(fft(x))/Nfft;
    fy1=ffy1(1:Nfft/2+1);
    fy1(2:end-1)=2*fy1(2:end-1);
    fx1=(0:Nfft/2)/Nfft*fs;
    figure
    subplot(3,1,1)
    plot(fx1,20*log(fy1));
    xlabel('频率/Hz');ylabel('FFT变换幅值/dB');title('信号的FFT变换');axis([0 80 -200 50]); 
    
    %Welch法功率谱估计
    [pxx1,w1]=pwelch(x,[],[],[],fs);
    subplot(3,1,2)
    plot(w1,10*log(pxx1))
    xlabel('频率/Hz');ylabel('功率密度/dB');title('信号的Welch法功率谱估计');axis([0 80 -200 50]); 
    %周期图法
    [pxx2,w2] = periodogram(x,[],[],fs);
    subplot(3,1,3)
    plot(w2,10*log(pxx2))
    xlabel(
    

    在这里插入图片描述

    前边谈到Welch法和周期图法的特点在上图中都有比较明显的体现,同时可以发现如果通过对比可以发现如果将FFT变换的频谱用指数形式表示,则最后的结果和周期图法得到的变换结果相一致,关于这一点相信大家也并不会意外。

    参数法功率图估计

    2. 信号的时频分析

    对于普通的FFT变换来说,其能够

    • 准确地反映信号所包含的频率分量,不能反映频率分量存在的时间段
    • 能反映调频信号的频率范围,不能反映频率随时间的变化规律
    • 能准确反映信号中某频率分量的总能量,缺不能显示频率分量强度强度随时间的变化功率

    简单地来说FFT是一种全局变换,在变换的过程中把整个信号的能量给平均了。我们知道FFT(x1+X2)=FFT(X1)+FFT(X2),傅里叶变换是满足线性性质的。如下图所示,对于信号1,2来说,这两个信号的FFT频谱一致,但信号而的余弦分量和信号1相比,以两倍的幅值出现了1半的时间。该谐波分量的总能量及平均幅值一样,故在FFT频谱上显示的幅值一致。

    对于平稳信号来说,FFT这么分析没有什么问题,但对于非平稳信号来说,如果想分辨1、2两种情况,仅仅用FFT变换就无法解决问题了,这时就需要借助这一节中需要将到的各类视频分析方法。
    在这里插入图片描述

    clear;clc;
    tic
    fs=100;
    T1=2;
    t1=(0:T1*fs-1)/fs;
    x1=cos(2*pi*15*t1)+0.2*rand(1,size(t1,2));
    T2=1;
    t2=(0:T2*fs-1)/fs;
    x2=2*cos(2*pi*15*t2)+0.2*rand(1,size(t2,2));
    x2=[x2,0.2*rand(1,size(t2,2))];
    Nfft=2^nextpow2(size(x1,2));
    %FFT变换
    ffy1=abs(fft(x1))/Nfft;
    fy1=ffy1(1:Nfft/2+1);
    fy1(2:end-1)=2*fy1(2:end-1);
    fx1=(0:Nfft/2)/Nfft*fs;
    
    ffy2=abs(fft(x2))/Nfft;
    fy2=ffy2(1:Nfft/2+1);
    fy2(2:end-1)=2*fy2(2:end-1);
    fx2=(0:Nfft/2)/Nfft*fs;
    %画图
    figure
    subplot(2,2,1)
    plot(t1,x1);
    xlabel('t/s');ylabel('幅值');title('时域信号1');axis([0,2,-3,3]);
    subplot(2,2,2)
    plot(t1,x2);
    xlabel('t/s');ylabel('幅值');title('时域信号2');axis([0,2,-3,3]);
    %FFT变换结果
    subplot(2,2,3)
    plot(fx1,fy1);
    xlabel('频率/Hz');ylabel('FFT变换幅值');title('信号1的FFT变换'); 
    subplot(2,2,4)
    plot(fx2,fy2);
    xlabel('频率/Hz');ylabel('FFT变换幅值');title('信号2的FFT变换'); 
    toc
    

    2.1 短时Fourier变换(Short Time Fourier Transform-STFT)

    STFT个人认为是其与传统的傅里叶变换之间并没有非常显著的区别。其本质上就是一个加窗的傅里叶变换,通过窗函数的移动来获取不同时间内的频率分量。傅里叶变换中拥有的性质基本上都可以直接移植到STFT上来。对于STFT来说,如果将窗口选得非常小,虽然时域上的分辨率非常高,但是由于时间太短,频域上的分辨率非常低。而如果窗口选得非常大的话,又变成了最原始的加窗傅里叶变换,时间上的分辨率非常低。时域和频域上的分辨率不可兼得。时域分辨率和频域分辨率的乘积是一个定值。
    在这里插入图片描述在这里插入图片描述在MATLAB中,官方给出的短时傅里叶变化你的函数为spectrogram,这里参考一篇博客和mathworks的官方文档,对该函数进行学习说明,该函数的用法为:
    [S,F,T,P]=spectrogram(x,window,noverlap,nfft,fs)
    在不给出左侧返回值的情况下,该函数直接生成使用短时傅里叶变换得到的频谱图。现对默认参数进行说明:

    • x-输入信号向量,默认被平分为8段,不能平分的话进行截断
    • window-若window为整数,则使用改长度的Hamming窗,若为向量,使用每段指定长度的窗函数
    • noverlap-各段之间重叠的采样点数,需要为小于window的整数
    • Nfft-进行离散傅里叶变换的点数
    • fs-采样频率,不指定的话默认为1hz
    • S-输入信号x的短时傅里叶变换,每一列包含一个短时局部时间的频率成分估计,时间沿列增加,频率沿行增加
    • F-??
    • T-频谱图计算的时刻点
    • P-能量谱密度???
      对于输入参数来说,nfft越大,频域的分辨率就越高,但时间的分辨率就越差,而noverlap影响时间轴分辨率,越大的话,时间分辨率越高
      下面使用一个具体的实例
    %计算并显示二次扫频信号的PSD图,扫频信号的频率开始于100Hz,在1s时经过200Hz
    T = 0:0.001:2;
    X = chirp(T,100,1,200,'q');
    figure()
    subplot(2,2,1);
    spectrogram(X,512,480,1024,1E3,'yaxis');
    title('512-Window-480-Overlap');
    subplot(2,2,2);
    spectrogram(X,512,120,1024,1E3,'yaxis');
    title('512-Window-120-Overlap');
    subplot(2,2,3);
    spectrogram(X,128,120,1024,1E3,'yaxis');
    title('128-Window-120-Overlap');
    subplot(2,2,4);
    spectrogram(X,128,10,1024,1E3,'yaxis');
    title('128-Window-10-Overlap');
    

    在这里插入图片描述从上面的对比中可以很明显的看出,当重叠较小的频谱在时域上看起来比较连续,但值得注意的是,由于这里频率变换比较快,且在指定时间频率成分较为单一,当窗函数较大时,指定窗函数内的频率成分较为多样,频率分辨率反而较低,需要根据具体的情况进行分辨。
    现使用实际数据来进行分析,取半个周期的数据

    clear;clc;
    %电机数据均为15Hz,内部负载为80load,采样率为5120Hz,取1个周期的数据
    %由于频率过高,进行降采样,频率为1280Hz,采样点为1280
    load('normal.mat');
    normal=data(1:4:5120*4,5);
    load('bowed.mat');
    bowed=data(1:4:5120*4,5);
    load('bearing.mat');
    bearing=data(1:4:5120*4,5);
    load('broken.mat');
    broken=data(1:4:5120*4,5);
    window=256;
    noverlap=200;
    nfft=2048;
    figure()
    subplot(2,2,1);
    spectrogram(normal,window,noverlap,nfft,5120,'yaxis');
    title('normal');
    subplot(2,2,2);
    spectrogram(broken,window,noverlap,nfft,5120,'yaxis');
    title('broken bar');
    subplot(2,2,3);
    spectrogram(bowed,window,noverlap,nfft,5120,'yaxis');
    title('bent rotor');
    subplot(2,2,4);
    spectrogram(bearing,window,noverlap,nfft,5120,'yaxis');
    title('bearing');
    

    在这里插入图片描述

    2.2 Gabor变换

    2.3 Wigner-Ville变换

    2.4 小波分析

    《连续小波变换、离散小波变换、二进小波变换、离散序列的小波变换、小波包》《有关小波的几个术语及常见的小波基介绍》

    简单总结

    • 连续小波变换
      就是,小波变换最常见到的公式和形式就是连续小波变换,公式如下:
      W f ( a , b ) = 1 ∣ α ∣ ∫ − ∞ ∞ f ( t ) Ψ ( t − b a ) d t W _ { f } ( a , b ) = \frac { 1 } { \sqrt { | \alpha | } } \int _ { - \infty } ^ { \infty } f ( t ) \Psi \left( \frac { t - b } { a } \right) \mathrm { d } t Wf(a,b)=α 1f(t)Ψ(atb)dt
      式子中a代表尺度(在某种意义上就是频率的概念),b是时间参数或者平移参数。不严谨的说 W f ( a , b ) W _ { f } ( a , b ) Wf(a,b)值得是对信号进行小波变换后当频率为a时间为b时的值,可以看出一维信号进过变换后成为了二维信号。
    • 离散小波变换
      离散小波变换DWT指的是将CWT中的尺度参数a和平移参数b离散化,而没有将信号和小波中的时间变量t离散化。
      在对尺度参数a和平移参数b离散的过程汇总,一般对尺度进行幂级数离散化,令 a = a 0 m a = a _ { 0 } ^ { m } a=a0m,对b进行均匀离散化,考虑到不同尺度下频率不同,因此不同尺度下参数b的离散间隔不同
    • 二进小波变换(Dyadic Wavelet Transform)
      在二进小波变换汇总,取特殊化,令 a 0 = 2 a_{0}=2 a0=2,然后保持平移参数b仍是连续的,则这类小波被称为二进小波变换。
    • 离散序列的小波变换
      上面说的三种小波变换的小波基都不是正交基,会带来一些麻烦,通过他们对信号进行变换后的信息是有冗余的,Mallat给出了一种在正交小波基上的信号分解算法,就是著名的Mallat算法。
      离散序列的小波变换就是基于Mallat算法,在分解的过程汇总,初始系数x(暂且这么称呼)与其第一层分解后的高频系数D1(细节部分Detail)的关系是x经过高通滤波器g滤波后再下采样,与低频系数A1(近似部分Approximate)的关系是x经过低通滤波器h滤波后再下采样;然后继续对低频系数A1进行第二层分解,低频系数A1与其第二层分解后的高频系数D2(细节部分Detail)的关系是A1经过高通滤波器g滤波后再下采样,与低频系数A1(近似部分Approximate)的关系是A1经过低通滤波器h滤波后再下采样;后面依次类推即可。由于一直在下采样,所以虽然滤波器系数g和h不变,但其滤波带宽一直在减半。初始系数是怎么来的呢?肯定是根据信号得到的,最简单最粗糙的办法就是对信号直接抽样。这是对连续信号进行正交小波分解,有了这些系数,再利用正交小波基,就可以表示出信号了,这类似于连续周期信号的傅里叶级数分解吧。
      在这里插入图片描述
    • 小波包
      无论是低频系数还是高频系统都进行同样的分解,然后选取一个最合适的分解路径。怎么评价分解是否是最优的呢?最自然的想法就是利益最大化或者是代价最小化,构建一个代价函数求一下看看如分解代价最小。代价函数有很多种,具体不说了。

    接下来使用MATLAB中的具体函数来对信号进行小波分解。

    这部分的内容较多,另外单独写了一篇文章MATLAB小波变换工具箱 Wavelet Toolbox 实际操作与训练,可以进行相关的查阅。

    不过小波变换虽然有诸多优点,但是小波变换始终受到小波基函数的影响,且一直没有非常合适的小波及选取理论,导致小波变换在信号处理的过程中无法实现对信号的自适应变换。

    2.5 Hilbert-Huang变换

    关于这部分内容可以参考一篇文章HHT变换讲义,个人觉得说的非常清楚,通俗易懂

    Hilbert-Huang变换,即HHT是一种用于处理不平衡,无规律数据的存在滋生调整型的数据处理方式。
    HHT主要内容包含两部分,第一部分为经验模态分解(Empirical Mode Decomposition,简称EMD),它是由Huang提出的;第二部分为Hilbert谱分析(Hilbert Spectrum Analysis,简称HSA)。简单说来,HHT处理非平稳信号的基本过程是:首先利用EMD方法将给定的信号分解为若干固有模态函数(以Intrinsic Mode Function或IMF表示,也称作本征模态函数),这些IMF是满足一定条件的分量;然后,对每一个IMF进行Hilbert变换,得到相应的Hilbert谱,即将每个IMF表示在联合的时频域中;最后,汇总所有IMF的Hilbert谱就会得到原始信号的Hilbert谱。
    在HHT中,IMF应该具有以下两条性质:

    • 极值总数与过零总数最多相差1
    • 上包络线和下包络线的均值必须是零

    模态分解的方法

    1.采用样条差值方法将x(t)的所有极大值点连线E1和极小值点连线E2,及平均包络线p(t)
    2.y(t)=x(t)-p(t),检查y是否满足IMF的性质,如果不满足的话反复操作,直到获得满足IMF性质的c(t)=y(t)作为x(t)的第1个IMF分量IMF1
    3.x-c得到全新的x,再重复1,2得到下一个IMF,直到满足规定要求后停止,x就能表示成为IMF和残差量q(t)的和
    x ( t ) = ∑ i = 1 m c i ( t ) + q ( t ) x ( t ) = \sum _ { i = 1 } ^ { m } c _ { i } ( t ) + q ( t ) x(t)=i=1mci(t)+q(t)

    Hilbert谱和Hilbert边际谱

    设序列x的Hilbert变换y为:
    y ( t ) = 1 π ∫ − ∞ ∞ x ( τ ) t − τ d τ = x ( t ) ∗ 1 π t y ( t ) = \frac { 1 } { \pi } \int _ { - \infty } ^ { \infty } \frac { x ( \tau ) } { t - \tau } d \tau = x ( t ) * \frac { 1 } { \pi t } y(t)=π1tτx(τ)dτ=x(t)πt1
    Hilbert变换的频率响应为
    H ( j ω ) = − sign ⁡ ( ω ) = { − j , ω ≥ 0 j , ω &lt; 0 H ( j \omega ) = - \operatorname { sign } ( \omega ) = \left\{ \begin{array} { c } { - j , \omega \geq 0 } \\ { j , \omega &lt; 0 } \end{array} \right. H(jω)=sign(ω)={j,ω0j,ω<0
    信号幅值没有变化,主要是改变了信号相位

    EMD在MATLAB中的应用

    早期的Matlab并没有HHT,EMD等内置函数,需要下载第三方的工具箱,或者自信编写相关的函数,而在R2018b的版本以后,MATLAB中已经有了相关的原生函数
    MATLAB中进行EMD分解的方法是:

    • [imf,residual,info] = emd(X,‘Interpolation’,‘pchip’);
      这种方法对信号进行EMD分解,并返回得到的固有模态函数IMF,残差量residual和分解过程中保存的信息info
      在这里插入图片描述
    • emd(X,‘Interpolation’,‘pchip’);
      如果没有左边参数,函数会直接返回直接画出原始信号,前三个IMF和残差量
      在这里插入图片描述
    • 参数说明:
      emd函数中可以对EMD分解过程的中的一些参数进行设置,如停止方式,可以选择最大迭代次数,误差收敛方式,最大IMF数量,等,具体可以参考Mathworks的官方文档。但Interpolation这个参数比较常用,需要注意, 如果信号x是光滑信号,应该选择spline,样条曲线,而如果信号x是非光滑的信号,则应该选择pchip,使用分段三次Hermite插值多项式方法进行连接。

    HHT变换

    对随机序列x(t)j进行希尔伯特变换之后得到y(t),由x,y可以得到一个解析信号z(t)
    Z ( t ) = X ( t ) + i Y ( t ) = a ( t ) e i θ ( τ ) Z ( t ) = X ( t ) + i Y ( t ) = a ( t ) e ^ { i \theta ( \tau ) } Z(t)=X(t)+iY(t)=a(t)eiθ(τ)
    其中 { a ( t ) = [ X 2 ( t ) + P 2 ( t ) ] 1 / 2 θ ( t ) = arctan ⁡ ( Y ( t ) X ( t ) ) \left\{ \begin{array} { l } { a ( t ) = \left[ X ^ { 2 } ( t ) + P ^ { 2 } ( t ) \right] ^ { 1 / 2 } } \\ { \theta ( t ) = \arctan \left( \frac { Y ( t ) } { X ( t ) } \right) } \end{array} \right. {a(t)=[X2(t)+P2(t)]1/2θ(t)=arctan(X(t)Y(t))
    在完成EMD后就提取了所有的IMF,对IMF进行Hilbert变换之后,就可以用下面的方式表示信号: X ( t ) = ∑ j = 1 n a j ( t ) exp ⁡ ( i ∫ ω j ( t ) d t ) X ( t ) = \sum _ { j = 1 } ^ { n } a _ { j } ( t ) \exp \left( i \int \omega _ { j } ( t ) d t \right) X(t)=j=1naj(t)exp(iωj(t)dt)
    上式中忽略了余项,因为余项的能量很小,一般可以忽略不计,由Hilbert变换得到的幅值和频率都是时间的函数,如果将幅值表示在频率时间的平面上,就可以得到Hilbert谱H(w,t),如果对H进行时间积分,就能够得到希尔伯特边际谱h(w)
    h ( w ) = ∫ 0 T H ( w , t ) d t h ( w ) = \int _ { 0 } ^ { T } H ( w , t ) d t h(w)=0TH(w,t)dt
    编辑普提供了对每个频率的总振幅的量测,表示了整个时间长度内累计的振幅,另外,作为Hilbert编辑普的附加结构,可以得到希尔伯特瞬时能量:
    I E ( t ) = ∫ w H 2 ( w , t ) d w I E ( t ) = \int _ { w } H ^ { 2 } ( w , t ) d w IE(t)=wH2(w,t)dw
    如果幅值的平方对时间积分,就可以得到Hilbe能量谱:
    E S ( w ) = ∫ 0 T H 2 ( w , t ) d t E S ( w ) = \int _ { 0 } ^ { T } H ^ { 2 } ( w , t ) d t ES(w)=0TH2(w,t)dt

    同emd函数一样,matlab也是在R2018b之后的版本中才有内置的hht函数,

    • hht(imf,fs);直接根据第一步获得的IMF画出Hilbert谱
      在这里插入图片描述
    • [hs,f,t] = hht(imf,fs)
      这里的返回值官方文档没有明确解释,对文献返回变量记性观察可以得到,t就是最简单的时间向量,f是将fs/2等分成100份后的一个101维的向量,而hs则是每个101*t维的矩阵,t是时间的变量的维度,同样也是IMF分量的函数

    这里如果不采用hht函数, 直接画图,而根据hht变换的原理,自行进行Hiltbert变换,会发现整体,过程较为繁琐下面列出部分代码

    z=hilbert(imf);
    instfreq = fs/(2*pi)*diff(unwrap(angle(z)));
    figure
    plot(t(2:end),instfreq)
    

    在这里插入图片描述
    上面只是画出了每个IMF在不同时间的频率分量,并没有表示强度,可以简单地尝试使用scatter来表示强度,代码及效果如下,但依然不是非常理想,还需要后续处理,这里暂时不做过多处理

    tm=repmat(t,1,9);
    for i=1:9
    scatter(t(2:end,i),instfreq(:,i),3,am(2:end,i),'filled')
    hold on
    end
    

    在这里插入图片描述

    分数阶Fourier变换

    其他一些问题说明

    是否转换为dB显示?

    这里转一小段文字来说明转换为分贝的作用。

    通常,振幅或功率谱以对数单位分贝(dB)的形式显示。该测量单位有助于查看宽动态范围,即可在存在较大信号分量时方便地查看小信号分量。分贝是比例单位,其计算方式如下。

    其中P是测量功率,Pr是参考功率。

    使用下列公式从振幅值计算分贝值。

    其中A是测量振幅,Ar是参考振幅。

    使用振幅或功率作为同一信号的振幅平方时,结果分贝水平是完全一致的。将分贝比乘以2,等同于将比例平方。因此,无论使用振幅或功率谱,都将得到相同的分贝水平和显示

    主要参考文献

    1. 《基于MATLAB的机械故障诊断技术案例教程》-张玲玲,肖静 高等教育出版社2016
      我们对Markdown编辑器进行了一些功能拓展与语法支持,除了标准的Markdown编辑器功能,我们增加了如下几点新功能,帮助你用它写博客:
    2. http://blog.sina.com.cn/s/blog_15183f5750102w13a.html
    3. https://www.penwatch.net/cms/fft_sampling/

    未完待续,持续更新中,欢迎大家指出相关问题,写到后面感觉前面突然一堆逻辑错误2333~

    留言回复不及时,大佬们烦劳发一下邮箱 zhao_lu_jie@163.com 哈,谢谢大家了!

    展开全文
  • 行业分类-电信-信号分析系统及方法.rar
  • 电力系统暂态信号的小波分析方法及其应用,MATLAB程序
  • 信号时域分析方法(有量纲特征值&无量纲特征值) 时域分析中无量纲特征值详解 时域特征值提取的Matlab实现 【频域】 信号频域分析方法 频域特征值提取的Matlab实现(频谱、功率谱、倒频谱) 频域特征值提取的...

    将不定期更新和补充
    以下链接均出自【括号先森】,为防止自己找不到,特加此汇总

    【时域】
    信号时域分析方法(有量纲特征值&无量纲特征值)
    时域分析中无量纲特征值详解
    时域特征值提取的Matlab实现

    【频域】
    信号频域分析方法
    频域特征值提取的Matlab实现(频谱、功率谱、倒频谱)
    频域特征值提取的Matlab实现(小波分析)

    以下转自知乎回答
    傅里叶变换与小波变换的关系

    【个人总结】
    《语音信号处理》书籍章节概览,可用于信号处理

    【Matlab代码实现——信号】
    matlab信号处理与滤波去噪
    matlab小波分析

    很有用

    【代码】(很有用)MATLAB仿真噪声信号、单多音信号、LFM、2ASK、2FSK、BPSK、16QAM
    【视频】(很有用)Matlab理论与实践视频(很有用):数字信号分析与实践

    视频教学

    【视频】b站算法工匠:matlab通信系统仿真(调制解调和信号同步)
    【视频】Matlabb站教学视频:Matlab系统仿真
    【视频】matlab官方视频:Matlab在信号处理中的应用
    【视频】b站视频:Matlab深度学习在深度学习的应用

    官方文档

    【官方网页】matlab官方:

    【官方包】

    【官方项目】

    matlab应用实践

    1. MATLAB仿真噪声信号、单多音信号、LFM、2ASK、2FSK、BPSK、16QAM
    2. I路Q路怎么画原始信号
    3. 如何用matlab仿真随机脉冲干扰
    4. 窄带干扰
    5. GPS信号和其干扰的Matlab仿真
    6. GPS信号和其干扰的Matlab仿真
    7. gps信号和干扰的matlab仿真代码
    8. matlab 产生伪随机序列
    9. 【网页文档】噪声调频、调相、调幅干扰实现
    10. 【网页文档】 线性调频信号MATLAB仿真
    11. 【网页文档】 QPSK调制与解调、I路和Q路信号
    12. 【网页文档】MATLAB不同时频信号处理方法介绍及效果对比
    13. 【网页文档】 MATLAB产生线性扫频波形
    14. 【网页文档】信号调制-3种方法:单一频率信号调制、OFDM信号调制、chirp信号调制
    展开全文
  • 信号的频域分析方法多种多样,这里针对较为常见的(频谱、能量谱、功率谱、倒频谱、小波分析)集中进行说明。这些方法的MATLAB代码实现参见文章频域特征值提取的MATLAB代码实现(频谱、功率谱、倒频谱) 在看这篇...

    信号的频域分析方法多种多样,这里针对较为常见的(频谱、能量谱、功率谱、倒频谱、小波分析)集中进行说明。这些方法的MATLAB代码实现参见文章频域特征值提取的MATLAB代码实现(频谱、功率谱、倒频谱)

    在看这篇文章之前可以参看之前的两篇,其中涉及一些时域特征值介绍和能量、功率信号区别的介绍:

    Mr.括号:信号的各种时域分析方法的理解

    Mr.括号:能量信号和功率信号的分别

    文章如要转载请私信与我联系,并注明来源知乎专栏与信号处理有关的那些东东作者Mr.括号。

     

    一、频谱

    一般来说,频谱分析指的是将信号做傅里叶变换从而进行分析。频谱分析是包括幅频谱和相频谱两张图的。不过最常用的是幅频谱。

    频谱分析是最常用和最重要也是最基础的频域分析方法。在这里不详细展开说了,不了解概念的可以参考这篇文章:

    Heinrich:如果看了这篇文章你还不懂傅里叶变换,那就过来掐死我吧

    时常会见到几个概念:DFT,DTFT,DFS,FFT,FT, FS等等。这些概念怎么区分?看这篇文章就能搞懂:一幅图弄清DFT与DTFT,DFS的关系 - BitArt - 博客园

    在这里概括一下(不感兴趣的话可以直接跳过):

    FT(Fourier Transformation):傅里叶变换。就是我们理论上学的概念,但是对于连续的信号无法在计算机上使用。其时域信号和频域信号都是连续的。

    DTFT(Discrete-time Fourier Transform):离散时间傅里叶变换。这里的“离散时间”指的是时域上式离散的,也就是计算机进行了采样。不过傅里叶变换后的结果依然是连续的。

    DFT(Discrete Fourier Transform):离散傅里叶变换。在DTFT之后,将傅里叶变换的结果也进行离散化,就是DFT。

    也就是说:FT时域连续、频域连续;DTFT时域离散、频域连续;DFT时域离散、频域离散。

    FFT(Fast Fourier Transformation):快速傅里叶变换。就是DFT的快速算法,一般工程应用时用的都是这种算法。

    FS(Fourier Series):傅里叶级数。是针对时域连续周期信号提出的,结果是离散的频域结果。

    DFS(Discrete Fourier Series):离散傅里叶级数。是针对时域离散周期信号提出的,DFS与DFT的本质是一样的。

    另外补充几点相关知识:

    • 在实际计算中通常使用快速傅里叶变换(FFT)。它是一种用来计算DFT(离散傅里叶变换)和IDFT(离散傅里叶反变换)的一种快速算法。

    • 随机信号是无法做傅里叶变换的(*这里要再补充)

     

    二、能量谱

    要理解能量谱和功率谱,首先要弄明白能量有限信号和功率有限信号(参看之前的文章能量信号和功率信号的分别)。

    能量谱也叫能量谱密度,能量谱密度描述了信号或时间序列的能量如何随频率分布。能量谱是原信号傅立叶变换的平方。

     

    三、功率谱

    功率谱是功率谱密度函数(PSD)的简称,它定义为单位频带内的信号功率。

    功率谱是针对功率信号来说的。功率谱的推导公式相对复杂,不过幸运的是维纳-辛钦定理证明了:一段信号的功率谱等于这段信号自相关函数的傅里叶变换。

    所以求功率谱就有了两种方法:1.(傅立叶变换的平方)/(区间长度);2.自相关函数的傅里叶变换。这两种方法分别叫做直接法和相关函数法。

    功率谱这里存在着一些问题,整理如下:

    1.功率谱密度的单位是什么,看有的写的是dB,还有的说是W/Hz。

    功率谱的单位是W/Hz,单位是dB时是做了对数处理(10logX)。取对数的目的是使那些振幅较低的成分相对高振幅成分得以拉高,以便观察掩盖在低幅噪声中的周期信号。

    2.求功率谱的两种方法有什么区别么?

    从原理上讲似乎没什么区别,从MATLAB仿真结果上来看,相关函数法对噪声的抑制效果更好,图线更平滑。(见频域特征值提取的MATLAB代码实现(频谱、功率谱、倒频谱)

    3.FFT和PSD都是表示的频谱特性,帮助我们找出峰值的位置,那么有了FFT为什么还要提出PSD?

    信号分为确定信号和随机信号,而确定信号又分为能量信号和功率信号,随机信号一定是功率信号。根据狄里赫利条件,能量信号可以直接进行傅里叶变换,而功率信号不行。对于无法做傅里叶变换的信号,只能走一步弯路,先求自相关,再做傅里叶。但是物理意义上就是功率谱了。不过总之得到了信号的频率特性。

    4.既然为什么随机信号的一次FFT没有意义却还能(傅立叶变换的平方)/(区间长度)得到功率谱?

    对随机信号直接做FFT的做法其实就是截断成能量信号进行处理,这种处理不符合随机信号定义,但之所以这样做,是做短时频域分析下作的近似处理。(这里希望能有大神能给出更好的解释)

    所以总结,频谱和能量谱(也叫能量谱密度)是傅里叶变换得到的复数结果和模平方的关系; 而功率谱(也就是功率谱密度)是针对随机信号分析提出的概念。

    更多讨论可以参看:随机信号傅里叶变换和功率谱密度图给出的信息有什么不同 - MATLAB中文论坛

     

    四、倒频谱

    倒频谱(Cepstrum)也叫倒谱、二次谱和对数功率谱等。倒频谱的工程型定义是:信号功率谱对数值进行傅立叶逆变换的结果。(信号→求功率谱→求对数→求傅里叶逆变换)

    为什么翻译作倒频谱呢?我个人的理解是,频谱(功率谱)反应的频率特征点横坐标是频率f(Hz),在倒频谱中对应的特征点的横坐标是时间t(s),而f与t互为倒数。从这里也可以看出,虽然倒频谱也叫“频谱”,其横坐标却并不是频率,而是时间。

    那么倒频谱有什么好处呢?

    “该分析方法方便提取、分析原频谱图上肉眼难以识别的周期性信号,能将原来频谱图上成族的边频带谱线简化为单根谱线,受传感器的测点位置及传输途径的影响小。”

    这都是啥意思?一条条解释:

    1.方便提取、分析原频谱图上肉眼难以识别的周期性信号

    我们知道,频谱分析就是为了提取原始信号中的周期性信号的,怎么频谱中的信号还会有周期性?这就又涉及到两个概念:调制和边频带。

    调制分为幅值调制和频率调制。下面以齿轮的幅值调制为例进行说明:齿轮的振动信号主要包括两部分,分别是齿轮啮合振动信号(高频)和齿轮轴的转频振动信号(低频),时域和频域曲线分别如下图所示:

    高频信号和低频信号时域波形

    高频信号和低频信号的频域波形

    调制就是高低频率信号的混合。幅值调制从数学上看,相当于两个信号在时域上相乘;而在频域上,相当于两个信号的卷积。调制后的信号在时域和频域上分别变为:

    调制后的时域信号

    调制后的频域信号

    我们发现,调制后的信号中,除原来的啮合频率分量外,增加了一对分量,它们是以高频信号特征频率为中心,对称分布于两侧,所以称为边频带。

    实际实验中齿轮啮合振动信号(高频)和齿轮轴的转频振动信号(低频)的特征频率可能是有多组的,其调制后的频域信号近似于一组频率间隔较大的脉冲函数和一组频率间隔较小的脉冲函数的卷积,从而在频谱上形成若干组围绕啮合频率及其倍频成分两侧的边频族,如下图:

    边频带的形成

    说了一大堆,终于回归到上边提到的问题:倒频谱“方便提取、分析原频谱图上肉眼难以识别的周期性信号”。这里指的周期性信号,就是重复出现的边频带。

    倒频谱能较好地检测出功率谱上的周期成分,通常在功率谱上无法对边频的总体水平作出定量估计,而倒频谱对边频成分具有“概括”能力,能较明显地显示出功率谱上的周期成分,将原来谱上成族的边频带谱线简化为单根谱线,便于观察,而齿轮发生故障时的振动频谱具有的边频带一般都具有等间隔(故障频率)的结构,利用倒频谱这个优点,可以检测出功率谱中难以辨识的周期性信号。

    2.受传感器的测点位置及传输途径的影响小

    这是倒频谱的第二个好处。对于布置在不同位置的传感器,由于传递路径不同,其功率谱也不相同。但在倒频谱上,由于信号源的振动效应和传递途径的效应分离开来,代表齿轮振动特征的倒频率分量几乎完全相同,只是低倒频率段存在由于传递函数差异而产生的影响。在进行倒频谱分析时,可以不必考虑信号测取时的衰减和标定系数所带来的影响。这一优点对于故障识别极为有用。

    关于倒频谱,文章 齿轮故障诊断常用信号分析处理方法 给出了具体了例子,方便理解。

     

    五、小波分析

    小波分析是一种时频域分析方法,该方法兼顾了信号在时域和频域的信息。知乎上有一篇文章对小波分析的理解进行了生动的讲解,强烈建议对小波分析概念不熟的同学先看一下。咚懂咚懂咚:能不能通俗的讲解下傅立叶分析和小波分析之间的关系?这篇文章中最后给出的小波变换的结果是这样的:

    连续小波变换

    看起来十分厉害,不过同时会发现两个问题:运算量很大;只有数值解,没有解析解。上述这种小波分析方法叫连续小波变换(continuous wavelet transform, CWT)。

    为了减少变换运算量,去除不必要的重复的系数,实际中使用的通常是离散小波变换(discrete wavelet transform, DWT)。

    这里的“离散”指的是什么呢?

    让我们先回到小波基波(也叫母小波)的表达式:

    其中s是尺度参数,表征频率;t是位移参数,表征时间。这部分在答友的连接里也提到了。再看上一张图,xy坐标分别是SCALE和TRANSLATION,也就是s和t,他们在连续小波变换中是连续的。

    所以,在离散小波变换中,“离散”的就是参数s和t。此时小波表达式写为:

    j和k都是整数,通常取s0=2,τ0=1。

    可以看出,随着j取值的递增,我们可以得到一串不同的小波(子小波,也叫女儿小波...)。这些子小波的尺度参数以2的j次方的形式增长。当使用这一系列的子小波,对一个连续函数进行离散分析时,我们所获得的是一组小波分析的系数,这个分析过程称为小波系列分解。

    上边说道,尺度参数表征的是频率,在子小波中尺度参数以2的倍数增长(即小波的“长度”被“拉长”了2倍),那么子小波对应能检测到的频率值也会以1/2的倍数缩小。母小波所对应的频谱位于频率谱的高端,具有最大的频率谱范围- 而其他的子小波的频率谱则依次向频谱图的低频端移动,同时它们所覆盖的频率谱范围也相应地递减。在理想的情况下,所有的滤波器应该首尾相接互相覆盖。

    不同尺度的子小波在小波频率谱上的覆盖

    是的,每个子小波就相当于一个滤波器,离散小波变换的过程就是逐级滤波的过程。

    具体流程是怎样的呢?

    用一句话描述就是:一组离散信号通过一系列的低通和高通滤波器,分别可以得到近似信号(用字母A表示)和细节信号(用字母D表示)。

    用一张图描述就是:

    LP为低通滤波器,HP为高通滤波器,B为带宽,2B为2倍带宽

    用一个例子来描述就是:

    对原始信号滤波:原始信号是在一个连续的低频正弦波信号(频率为0.5)上随机叠加了两个高频(频率为 10)高振幅的正弦脉冲,这里使用了 dB5(第五级Daubechies小波)作为去除噪音操作的母小波。

    原始信号

    一阶小波分解的结果为:

    一阶小波分解的近似信号(低通结果)

    一阶小波分解的细节信号(高通结果)

    二阶小波分解的结果为(即对A1信号做分解):

    二阶小波分解的近似信号(低通结果)

    二阶小波分解的细节信号(高通结果)

    三阶小波分解的结果为(即对A2信号做分解):

    三阶小波分解的近似信号(低通结果)

    三阶小波分解的细节信号(高通结果)

    四阶小波分解的结果为(即对A3信号做分解):

    四阶小波分解的近似信号(低通结果)

    四阶小波分解的细节信号(高通结果)

    至此我们已经能够得到较好的滤波结果了(即A5,不过脉冲信号也被滤掉了,用设置门限的方法可以保留住该信号,这里不做展开)。可以看到原式信号被逐级的,无遗漏地进行了高、低通滤波,且越接近低频分段越细,几乎想要哪个频段的特征都能得到,因而这个方法有个霸气的名字,叫滤波器银行。

    小波变换大致讲完了,那么它有哪些特点,可以用来做什么呢?

    首先当然是滤除噪声。

    其次,由于离散小波变换是可逆的,所以还可以用来做图像、信号的无损压缩。

    另外,可以检测信号的非连续性。奇异信号会在细节信号D中展现。

    除此之外,还可以进行图像边界检测等工作。

     

    (时)频域的一些常见的分析处理方法就介绍到这里了。这里能讲到的只是一些表面的概念和一些个人的理解。文中的内容还会不断修正和完善,对于不严谨以及谬误指出,还望多多指教。

    欢迎关注我的公众号“括号的城堡”,微信号为“khscience”,会有更多有趣的东西分享。

     

    参考:

    Heinrich:如果看了这篇文章你还不懂傅里叶变换,那就过来掐死我吧

    随机信号傅里叶变换和功率谱密度图给出的信息有什么不同 - MATLAB中文论坛

    傅里叶变换和逆变换公式的我理解意义 - CSDN博客

    齿轮的振动机理

    齿轮故障诊断常用信号分析处理方法

    Cepstrum - Wikipedia

    Wavelet - Wikipedia

    周宇峰, 程景全. 小波变换及其应用[J]. 物理, 2008, 37(1):24-32.

    展开全文
  • 行业分类-电信-一种信号分析装置及其配置方法.rar
  • 基于数字信号处理器的12通道心电图机的设计及分析方法的研究
  • 章立军博士学位论文《信号的数学形态学分析方法及其应用研究》,主要研究了如何将数学形态学应用于故障诊断领域。
  • 参数化分析方法有ARMA模型、分形维数计算用于判别微弱信号存在性的混沌分析法等;非参数化方法即谱分析方法,包括传统的DFT谱方法在此基础上的提高频率分辨能力(如AR功率谱)和时频联合分析(小波分解、Wigner...
  • MATLAB信号处理地震波当中的应用。大家注意分享哦,呵呵,有好的我会继续上传的。
  • 信号进行相位噪声指标测量是现在工作中经常遇到的事情,本文首先从信号相位噪声的定义入手,重点介绍使用信号分析仪进行相位噪声测量的方法及注意事项。  1、相位噪声是什么?  在频域内,一个理想正弦波信号...
  • 介绍了一种处理非线性、非平稳信号的新方法――HHT的原理特点,并将其应用于不同思维作业脑电信号分析。实验结果表明,不同思维作业脑电信号经HHT后的HH谱和Hilbert边际谱都差异显著,证明HHT方法对脑电信号处理的...
  • 针对强烈的非线性、非Gauss性和非平稳性的旋转机械振动信号,提出一类基于小波包分析技术和传统的时频分析技术相结合的非平稳信号的短时快变过程时频分析方法,利用小波变换来实现变尺度的时频原子(小波族函数)信号...
  • 岸基GPS海面反射信号提取方法及特性分析
  • 采用理论分析、MATLAB模拟仿真和现场试验相结合的方法,对比分析了谐波权重测频法与汉宁窗傅里叶测频法、直接傅里叶测频法和硬件边沿测频法等方法不同频率的测试精度和误差。结果表明:在有限采样样本高实时性的...
  • 爆破振动信号分析技术是爆破危害...根据IMF分量经过Hilbert变换后得到的Hilbert能量谱,从时间-频率-能量上研究振动信号不同频率产生的影响,验证了HHT分析方法在非线性非稳态的爆破信号中的高效性和良好的自适应性。
  • 行业分类-电信-基于振动信号频谱分析的GIS故障诊断系统及方法.rar
  • 对135发动机纯压缩缸内压力两种实测信号和由KIVA2计算的模拟信号的三者频谱对比分析,得到了通道效应干扰在发动机缸内压力频谱中的事实位置形状,并通过有、无通道效应的两种实测燃烧压力时域信号频谱对比分析得到...
  • 岸基GPS海面反射信号提取方法及特性分析.pdf
  • 行业分类-电信-具有莫尔斯信号分析及转换功能的移动通信终端及方法.rar
  • 信号及频谱分析

    千次阅读 2020-04-26 13:15:07
    2、从分析域上分:时域信号、频域信号 3、从信号波形分:连续时间信号、离散时间信号 4、连续时间信号又可以分为 :动态信号信号的幅值、相位、周期等特征参数随时间的变化而变化的信号)、静态信号 动态信号...

    一、信号的分类

    1、从信号描述上分:确定性信号(可以用明确数学关系来描述)、非确定性信号

    2、从分析域上分:时域信号、频域信号

    3、从信号波形分:连续时间信号、离散时间信号

    4、连续时间信号又可以分为 :动态信号(信号的幅值、相位、周期等特征参数随时间的变化而变化的信号)、静态信号

                    动态信号又可以分为:确定性信号非确定性信号

                                               确定性信号可以分为:周期信号(经过一定时间可以重复出现的信号)、非周期信号

                                               非确定性信号可以分为:平稳随机信号非平稳随机信号

                                                            周期信号分为:简单周期信号(单一频率)、复杂周期信号(多个频率)

                                                            非周期信号分为:准周期信号(由多个周期信号合成,但各信号频率不成公倍数)、瞬态信号

    二、时域分析与频域分析的概念

    1、时域描述:信号用幅值随时间的变化来表示,通常称为时域分析(波形分析)。

    时域分析只能反映信号的幅值随时间的变化情况,除单频率分量的简谐波外,很难明确揭示信号的频率组成和各频率分量大小。

    2、频域分析:以频率为横坐标描述信号的频率结构和频率成分的幅值、相位关系。

    3、频域分析的目的:为了研究信号的频率结构和各频率成分的幅值、相位关系,应对信号进行频谱分析,把信号的时域描述通过适当方法变成信号的频域描述,以频率为独立变量来表示信号。

    如图所示,对复杂信号进行频域分析,可以活得各频率分量的幅值和相位关系。

    三、频域分析实例

    1、方波信号

    2、用傅立叶级数展开该周期方波信号

    w0处的幅值为       4*A/(Π),相位为0,

    3*w0处的幅值为   4*A/(3*Π),相位为0,

    ……

    7*w0处的幅值为   4*A/(7*Π),相位为0,

    根据傅里叶级数展开表达式可得不同频率下的幅值和相位角,并绘制出幅频图(右上)和相频图(中下)。

    3、方波信号波德图


    4、周期信号的频谱分析

    定义:信号频谱分析是采用傅立叶变换将时域信号x(t)变换为频域信号X(f),从而帮助人们从另一个角度来了解信号的特征。

    对于任何一个周期为T、且定义在区间(- T/2, T/2)内的周期信号f(t),都可以用上述区间内的三角傅立叶级数表示:

    式中的系数为:

    将方波信号按照傅里叶级数展开可得,

    根据计算结果,可画出方波信号的幅频图和相频图

    工程上习惯将计算结果用图形方式表示,以ωn为横坐标,幅值、相位角为纵坐标画图,则称为幅值-相位谱。

    傅里叶展开式的另一种形式

    四、相位的概念

    1、什么是相位

    相位只要是针对交流电信号来讲的,是描述交流电信号的特征参数之一。相位是反映交流电任何时刻的状态的物理量

    例如,对于正弦交流电信号I=A*sin(2*Π*f*t)来说,y的大小和方向是随时间变化的

    I是交流电流的瞬时值,A是交流电流的最大值,f是交流电的频率,t是时间。随着时间t的变化,交流电流可以从零变到最大值,从最大值变到零,又从零变到负的最大值,从负的最大值变到零。在三角函数中2πft相当于角度,它反映了交流电任何时刻所处的状态,是在增大还是在减小,是正的还是负的等等。因此把2πft叫做相位,或者叫做相。

    如果t等于零的时候,I并不等于零,公式应该改成I=Isin(2πft+ψ)。那么2πft+ψ叫做相位,ψ叫做初相位,或者叫做初相。

    相位(phase)是对于一个波形来说,能反映它在特定的时刻的位置:一种它是否在波峰波谷或它们之间的某点的标度。是描述讯号波形变化的度量,通常以度(角度)作为单位,也称作相角。

    2、相位的超前与滞后

    如图所示的信号分别为y=sin(2*pi*t);y1=sin(2*pi*t+pi/2);

    y1的波峰滞后于y的波峰,我们称y1的相位滞后于y的相位,滞后的相位差为phi=(t2-t1)*360/T(T为周期)。也可以说y1超前于y,超前的相位差为360-phi。

    五、MATLAB之傅里叶变换函数FFT

    1、变换流程

    初始信号为y=10*sin(2*pi*1*t+phi),式中1为固有频率,phi为初相位。

    确定采集信号的序列长度:

    对连续信号的采样,确定采样频率fs=1000

    确定数据长度 N=512。

    横坐标t的从时域到频域的变换:

    其变换公式为:式中,采样时间为  ts,数据长度为 N,任意数字频率为 k,

    幅值y从时域到频域的变换:

    Y=fft(y,N); %对 y 进行傅里叶变换 ,N为数据长度

    幅频图Y的绘制

    plot(fk,Y)

    相频图phi的绘制

    phi=angle(Y(1:k+1))*180/pi %求得 k+1个幅值点对应的相位角

    plot(fk,phi)

    2、matlab代码

    clc
    clear
    A=10;
    fw=1;%固有频率
    phi=pi/3;%相位
    fs=1000; %采样频率
    t=0:1/fs:10*pi;%采样时间范围t   采样时间为ts=1/1000=1/step
    y=A*sin(2*pi*fw*t+phi);%正弦函数 y
    
    f=fs*(0:256)/512;%(0:256)表示257个谱线,512为数据长度 采样时间为ts=1/1000=1/step
    %采样时间范围>>对应工程频率(采样时间转换为频率) 
    subplot(3,1,1);%三行一列第一幅图
    plot(t,y);%绘制图形
    xlabel('t/s','fontsize',13);
    ylabel('y','fontsize',13);
    title('正弦函数曲线','fontsize',13);%显示标题
    Y=fft(y,512);%对 y 进行傅里叶变换 取512个数据
    subplot(3,1,2);%三行一列第二幅图
    plot(f,abs(Y(1:257)));%绘制图形 f为频域,Y(1:257)为幅频谱的257个点
    xlabel('f/Hz','fontsize',13);%横坐标显示 f/Hz,字号 13 
    ylabel('幅值','fontsize',13);%纵坐标显示幅值,字号 13 
    title('幅频特性曲线','fontsize',13);%显示标题 
    [value,index]=max(abs(Y));%将 abs(Y)最大值点的横
    
    subplot(3,1,3);%三行一列第三幅图
    plot(f,angle(Y(1:257))*180/pi);%绘制图形
    xlabel('f/Hz','fontsize',13);%横坐标显示 f/Hz,字号 13
    ylabel('相位/°','fontsize',13);%纵坐标显示相位/°,字号 13
    title('相位特性曲线','fontsize',13);%显示标题
    

     

    展开全文

空空如也

空空如也

1 2 3 4 5 ... 20
收藏数 139,541
精华内容 55,816
关键字:

不同信号及其分析方法