精华内容
下载资源
问答
  • 研究结果表明,在共轭傅里叶变换光谱重构中,引入Nuttall能够有效减小频谱泄露对光谱分辨率的影响,使光谱失真得到有效抑制。与未引入Nuttall相比,该方法获得的复原光谱与原始光谱更加吻合。
  • 短时傅里叶变换(Short-timeFouriertransform,STFT)是一种常见的时频分析工具,窗函数的选择是其分析的前提。通过Matlab仿真讨论了在软件解调中不同分析窗函数及窗口长度对解调性能的影响,针对软件解调中解调精度和...
  • [TFR,T,F]=TFRSTFT(X,T,N,H,TRACE)  X : 信号. ... H : 选择的平滑窗函数 (默认值: Hamming(N/4)).  TRACE : 如果不为零,跟踪计算过程 (默认为 : 0). TFR : 得到的时频谱值(复数)  F : 归一化...

    [TFR,T,F]=TFRSTFT(X,T,N,H,TRACE)
      X : 信号.
      T : 时间序列 (默认值 :1:length(X)).
      N : 频率点数 (默认值: length(X)).
      H : 选择的平滑窗函数 (默认值: Hamming(N/4)).
      TRACE : 如果不为零,跟踪计算过程 (默认为 : 0).
    TFR : 得到的时频谱值(复数)
      F : 归一化频率(-0.5 — +0.5),实际频率FF=fs*F;

    clc
    clear all
    fs=1024;
    N=2048;
    t1=(0:N-1)/fs
    T=(0:N/2-1)/fs
    a=sin(2*pi*5*t1);
    subplot(311)
    plot(t1,a);
    [tfr,t,f] = tfrstft(a');%a只能有一列
    freq=(0:N/2)*fs/N;
    subplot(312)
    plot(freq,abs(tfr(1:N/2+1)));
    k=fft(a)
    subplot(313)
    plot(freq,abs(k(1:N/2+1)));
    

    在这里插入图片描述

    展开全文
  • 窗函数-减少傅里叶变换泄漏

    万次阅读 2017-01-12 10:53:44
    数字信号处理的主要数学工具是傅里叶变换.而傅里叶变换是研究整个时间域和频率域的关系。不过,当运用计算机实现工程测试信号处理时,不可能对无限长的信号进行测量和运算,而是取其有限的时间片段进行分析。做法是...

    数字信号处理的主要数学工具是傅里叶变换.而傅里叶变换是研究整个时间域频率域的关系。

    快速傅里叶变换假定了时间信号是周期无限的。但在分析时,我们往往只截取其中的一部分,因此需要加窗以减小泄露。窗函数可以加在时域,也可以加在频域上,但在时域上加窗更为普遍。

    不过,当运用计算机实现工程测试信号处理时,不可能对无限长的信号进行测量和运算,而是取其有限的时间片段进行分析。做法是从信号中截取一个时间片段,然后用截取的信号时间片段进行周期延拓处理,得到虚拟的无限长的信号,然后就可以对信号进行傅里叶变换、相关分析等数学处理。无限长的信号被截断以后,其频谱发生了畸变,原来集中在f(0)处的能量被分散到两个较宽的频带中去了(这种现象称之为频谱能量泄漏)。

    几种常用的窗函数的比较

    名称

    特点

    应用

    矩形窗

    Rectangle

    矩形窗使用最多,习惯上不加窗就是使信号通过了矩形窗。这种窗的优点是主瓣比较集中,缺点是旁瓣较高,并有负旁瓣,导致变换中带进了高频干扰和泄漏,甚至出现负谱现象。频率识别精度最高,幅值识别精度最低,所以矩形窗不是一个理想的窗。

    如果仅要求精确读出主瓣频率,而不考虑幅值精度,则可选用矩形窗,例如测量物体的自振频率等,也可以用在阶次分析中。

    汉宁窗

    Hanning

    又称升余弦窗。主瓣加宽并降低,旁瓣则显著减小,从减小泄漏观点出发,汉宁窗优于矩形窗.但汉宁窗主瓣加宽,相当于分析带宽加宽,频率分辨力下降。它与矩形窗相比,泄漏、波动都减小了,并且选择性也提高。

    是很有用的窗函数。如果测试信号有多个频率分量,频谱表现的十分复杂,且测试的目的更多关注频率点而非能量的大小,需要选择汉宁窗。如果被测信号是随机或者未知的,选择汉宁窗。

    海明窗

    (汉明窗)

    Hamming

    与汉宁窗都是余弦窗,又称改进的升余弦窗,只是加权系数不同,使旁瓣达到更小。但其旁瓣衰减速度比汉宁窗衰减速度慢。

    与汉明窗类似,也是很有用的窗函数。

    平顶窗

    Flap Top

    平顶窗在频域时的表现就象它的名称一样有非常小的通带波动。

    由于在幅度上有较小的误差,所以这个窗可以用在校准上。

    凯塞窗

    Kaiser

    定义了一组可调的由零阶贝塞尔Bessel 函数构成的窗函数,通过调整参数β可以在主瓣宽度和旁瓣衰减之间自由选择它们的比重。对于某一长度的Kaiser 窗,给定β,则旁瓣高度也就固定了。

     

    布莱克曼窗

    Blackman

    二阶升余弦窗,主瓣宽,旁瓣比较低,但等效噪声带宽比汉宁窗要大一点,波动却小一点。频率识别精度最低,但幅值识别精度最高,有更好的选择性。

    常用来检测两个频率相近幅度不同的信号。

    高斯窗

    Gaussian

    是一种指数窗。主瓣较宽,故而频率分辨力低;无负的旁瓣,第一旁瓣衰减达一55dB。常被用来截短一些非周期信号,如指数衰减信号等。

    对于随时间按指数衰减的函数,可采用指数窗来提高信噪比。

    三角窗

    (费杰窗)

    Fejer

    是幂窗的一次方形式。与矩形窗比较,主瓣宽约等于矩形窗的两倍,但旁瓣小,而且无负旁瓣。

    如果分析窄带信号,且有较强的干扰噪声,则应选用旁瓣幅度小的窗函数,如汉宁窗、三角窗等;

    切比雪夫窗(Chebyshev)

    在给定旁瓣高度下,Chebyshev窗的主瓣宽度最小,具有等波动性,也就是说,其所有的旁瓣都具有相等的高度。

     

    窗函数是频谱分析中一个重要的部分,窗函数修正了由于信号的非周期性并减小了频谱中由于泄露而带来的测量不准确性

    快速傅里叶变换假定了时间信号是周期无限的。但在分析时,我们往往只截取其中的一部分,因此需要加窗以减小泄露。窗函数可以加在时域,也可以加在频域上,但在时域上加窗更为普遍。截断效应带来了泄漏,窗函数是为了减小这个截断效应,其设计成一组加权系数。例如,一个窗函数可以定义为:
    w(t)=g(t)      -T/2<t<T/2
    w(t)=0         其他

    g(t)是窗函数,T是窗函数的时间.

    待分析的数据x(t)则表示为:

    x(t)=w(t)*x(t)'

    x(t)'表示原始信号,x(t)表示待分析信号。

     

    加窗在时域上表现的是点乘,因此在频域上则表现为卷积。卷积可以被看成是一个平滑的过程。这个平滑过程可以被看出是由一组具有特定函数形状的滤波器,因此,原始信号中在某一频率点上的能量会结合滤波器的形状表现出来,从而减小泄漏。基于这个原理,人们通常在时域上直接加窗。

    大多数的信号分析仪一般使用矩形窗(rectangular),汉宁(hann),flattop和其他的一些窗函数。

    矩形窗函数:
    w(k)=1

    汉宁窗:     
    w(k)=0.5*(1-cos(2*pi*k/(N-1)))    0<=k<=N-1

    由于加窗计算中衰减了原始信号的部分能量,因此对于最后的结果还需要加上修正系数。在线性谱分析中,一般使用幅度系数(amplitudecorrection),在功率谱中,一般使用能量系数(energy correction)。(这段不清楚在实际中如何用)

     

    matlab中提供了很多窗函数,如下

    image

    还提供了显示窗函数的GUI工具,如wvtool可以显示用来显示窗的形状和频域图形,wintool可以打开窗设计和分析工具,如运行

    wvtool(hamming(64),hann(64),gausswin(64))

    可以对比汉明窗、汉宁窗和高斯窗

    image

     

    简单测试一下加窗的效果如下

    image

    image

    可以看到加窗后,频谱泄露确实减少了,但同时信号能量也减小了,这也许就是所说的要使用能量系数吧,如下,这样一来,对比就更明显了,加窗可以有效的减少频谱泄露。

    image

     

    测试代码如下

    %% 窗函数测试
    function main
    clc
    close all

    Ts = 0.001;
    Fs = 1/Ts;
    %% 原始信号
    t = 0:Ts:pi/2;
    yt = sin(2*pi*5*t) + sin(2*pi*10*t) + sin(2*pi*15*t);
    [Yf, f] = Spectrum_Calc(yt, Fs);

    figure
    subplot(211)
    plot(t, yt)
    xlabel('t')
    ylabel('y')
    title('原始信号')
    subplot(212)
    plot(f, Yf)
    xlabel('f')
    ylabel('|Yf|')
    xlim([0 100])
    ylim([0 1])
    title('原始信号频谱')

    %% 加窗信号
    win = hann(length(t));
    yt1 = yt.*win';
    [Yf1, f1] = Spectrum_Calc(yt1, Fs);

    figure
    subplot(211)
    plot(t, yt1)
    xlabel('t')
    ylabel('y')
    title('加窗信号')
    subplot(212)
    plot(f1, 2*Yf1) % 2表示能量系数
    xlabel('f')
    ylabel('|Yf|')
    xlim([0 100])
    ylim([0 1])
    title('加窗信号频谱')
    end

    %% 求取频谱
    function [Yf, f] = Spectrum_Calc(yt, Fs)
    L = length(yt);

    NFFT = 2^nextpow2(L);
    Yf = fft(yt,NFFT)/L;

    Yf = 2*abs(Yf(1:NFFT/2+1));
    f = Fs/2*linspace(0,1,NFFT/2+1);
    end


    展开全文
  • 加窗傅里叶变换的演示 matlab程序 分别对加方窗和海明的信号做傅里叶变换 函数可以改变窗口的大小 对理解傅里叶变换和频谱非常有帮助
  • 离散傅里叶变换公式: X(k)=∑n=0N−1xne−j2πknN X(k)=\sum_{n=0}^{N-1}x_ne^{-j\frac{2\pi kn}{N}} X(k)=n=0∑N−1​xn​e−jN2πkn​ 下文中将 e−j2πknNe^{-j\frac{2\pi kn}{N}}e−jN2πkn​ 记为WNknW_N^{kn}...

    离散傅里叶与卷积

    离散傅里叶变换公式:
    X(k)=n=0N1xnej2πknN X(k)=\sum_{n=0}^{N-1}x_ne^{-j\frac{2\pi kn}{N}}
    下文中将 ej2πknNe^{-j\frac{2\pi kn}{N}} 记为WNknW_N^{kn} 也就是DFT变换公式下文记为
    X(k)=n=0N1xnWNkn X(k)=\sum_{n=0}^{N-1}x_nW^{kn}_{N}
    那么离散傅里叶逆变换IDFT表示为:
    Xn=k=0N1X(k)WNkn Xn=\sum_{k=0}^{N-1}X(k)W^{-kn}_{N}
    离散数字的卷积公式:
    xn=m=0M1u(nm)v(m) ,vM()uNnmu() x_n=\sum_{m=0}^{M-1} u(n-m)v(m) \ ,v为M个数据序列(卷积核),u为N个数据序列 \\ 当n-m小于零时,u补零(等于零)

    离散数字的卷积与傅里叶变换的关系:
    卷积后的DFT为:
    X(k)=n=0N1xnWNkn=n=0N1WNkn(m=0M1u(nm)v(m))=n=0N1(m=0M1u(nm)v(m)WNkn) \begin{aligned} X(k) &= \sum_{n=0}^{N-1}x_nW^{kn}_{N} \hspace{8cm} \\ &= \sum_{n=0}^{N-1}W^{kn}_{N}\left( \sum_{m=0}^{M-1} u(n-m)v(m) \right) \\ &= \sum_{n=0}^{N-1}\left( \sum_{m=0}^{M-1} u(n-m)v(m)W^{kn}_{N} \right) \\ \end{aligned}
    上面最后化简公式从m维度先求和 跟从n维度先求和 结果是一样的 也就是说
    n=0N1(m=0M1u(nm)v(m)WNkn)=m=0M1(n=0N1u(nm)v(m)WNkn) \begin{aligned} \sum_{n=0}^{N-1}\left( \sum_{m=0}^{M-1} u(n-m)v(m)W^{kn}_{N} \right) \hspace{4cm} \\ =\sum_{m=0}^{M-1}\left( \sum_{n=0}^{N-1} u(n-m)v(m)W^{kn}_{N} \right) \hspace{4cm} \end{aligned}
    继续化简
    所以有卷积后的DFT为:
    X(k)=m=0M1(n=0N1u(nm)v(m)WNkn)=m=0M1(n=0N1u(nm)v(m)WNkn)=m=0M1(n=0N1u(nm)v(m)WNk(nm)WNkm)=m=0M1v(m)WNkm(n=0N1u(nm)WNk(nm)) \begin{aligned} X(k) &= \sum_{m=0}^{M-1}\left( \sum_{n=0}^{N-1} u(n-m)v(m)W^{kn}_{N} \right) \hspace{4cm} \\ &= \sum_{m=0}^{M-1}\left( \sum_{n=0}^{N-1} u(n-m)v(m)W^{kn}_{N} \right) \\ &= \sum_{m=0}^{M-1}\left( \sum_{n=0}^{N-1} u(n-m)v(m)W^{k(n-m)}_{N}W^{km}_{N} \right) \\ &= \sum_{m=0}^{M-1}v(m)W^{km}_{N}\left( \sum_{n=0}^{N-1} u(n-m)W^{k(n-m)}_{N} \right) \end{aligned}
    化简结果跟X(k)=n=0N1xnWNknX(k)=\sum_{n=0}^{N-1}x_nW^{kn}_{N} 是有相似的形式的,
    如果化简的后半部分是nm=0N1u(nm)WNk(nm)\sum_{n-m=0}^{N-1} u(n-m)W^{k(n-m)}_{N},就可以认为离散数据卷积后
    的DFT等于卷积前两函数的DFT相乘了,但是实际上不是啊
    那么,n=0N1u(nm)WNk(nm)\sum_{n=0}^{N-1} u(n-m)W^{k(n-m)}_{N}nm=0N1u(nm)WNk(nm)\sum_{n-m=0}^{N-1} u(n-m)W^{k(n-m)}_{N} 有没有什么关系呢?
    当m=0时,它们是相等的,
    当m>0 ,n-m<0时 u(n-m)恒等于0 ,且令 n>N-1 u(n-m)也等于0
    有n-m<0 和 n>N-1 时 WNk(nm)×0=0W_{N}^{k(n-m)}\times 0=0
    则下面的等式是成立的
    n=0N1u(nm)WNk(nm)=n=mm+N1u(nm)WNk(nm)=U(k) \begin{aligned} \sum_{n=0}^{N-1} u(n-m)W^{k(n-m)}_{N} &= \sum_{n=m}^{m+N-1} u(n-m)W^{k(n-m)}_{N} \\ &=U'(k) \end{aligned}
    U(k)U'(k)是从原来u数据序列抽取第0,1,2…N-1-m; N-m个数,再补m个零组成一个新的N个数据u’的DFT结果,
    显然U’(k)与原来u数据序列DFT结果U(k)是不相等的,但当N远远大于m时,U’(k)与U(k)是差不多相等的;因为从一组N个数据中抽取大部分数,补零后重新组成新的N个数据,两者的频普是差不多相等的
    ,也可以说DFT是差不多相等的
    所以对上面的公式进一步化简

    X(k)=m=0M1v(m)WNkm(n=0N1u(nm)WNk(nm))=m=0M1v(m)WNkm(n=mm+N1u(nm)WNk(nm))m=0M1v(m)WNkmU(k)U(k)m=0M1v(m)WNkmU(k)V(k) \begin{aligned} X(k) &= \sum_{m=0}^{M-1}v(m)W^{km}_{N}\left( \sum_{n=0}^{N-1} u(n-m)W^{k(n-m)}_{N} \right) \\ &= \sum_{m=0}^{M-1}v(m)W^{km}_{N}\left( \sum_{n=m}^{m+N-1} u(n-m)W^{k(n-m)}_{N} \right) \\ & \approx \sum_{m=0}^{M-1}v(m)W^{km}_{N}U(k) \\ & \approx U(k)\sum_{m=0}^{M-1}v(m)W^{km}_{N} \\ & \approx U(k)V(k) \end{aligned}
    这里还有一个问题,就是v没有与u相当的频率数量,要让上面的约等式成立,
    也必须补零扩展v
    v补零不会对上面说到的N远远大于m条件冲突,因为当m>M-1时,v(m)为零,相当于m=0M1v(m)WNkmU(k)\sum_{m=0}^{M-1}v(m)W^{km}_{N}U(k)多加了若干个零。(但是v补零会使V(k)引入一些无关的频率,从而使结果误差变大)
    X(k)m=0M1v(m)WNkmU(k)m=0N1v(m)WNkmU(k)U(k)m=0N1v(m)WNkmU(k)V(k) \begin{aligned} X(k) & \approx \sum_{m=0}^{M-1}v(m)W^{km}_{N}U(k) \\ & \approx \sum_{m=0}^{N-1}v(m)W^{km}_{N}U(k) \\ & \approx U(k)\sum_{m=0}^{N-1}v(m)W^{km}_{N} \\ & \approx U(k)V(k) \end{aligned}

    结论:对离散数据卷积,相当于其频域与卷积核的域相乘,也就是起到滤波的效果


    离散数字的乘积与傅里叶变换后卷积的关系:

    对两组数据 u,v DFT后的U(k) V(k)进行卷积,再进行傅里叶逆变换:
    xn=k=0N1WNkn(m=0M1U(km)V(m))=m=0M1V(m)WNnm(k=0N1U(km)WNn(km))=m=0M1V(m)WNnm(k=mm+N1U(km)WNn(km))m=0N1v(m)WNnmununm=0N1v(m)WNnmunvn...DFT(unvn)=n=0N1unvnWNknm=0N1U(km)V(m) \begin{aligned} x_n &= \sum_{k=0}^{N-1}W^{-kn}_{N}\left( \sum_{m=0}^{M-1} U(k-m)V(m) \right) \\ &= \sum_{m=0}^{M-1}V(m)W^{-nm}_{N}\left( \sum_{k=0}^{N-1} U(k-m)W^{-n(k-m)}_{N} \right) \\ &= \sum_{m=0}^{M-1}V(m)W^{-nm}_{N}\left( \sum_{k=m}^{m+N-1} U(k-m)W^{-n(k-m)}_{N} \right) \\ & \approx \sum_{m=0}^{N-1}v(m)W^{-nm}_{N}u_n \\ & \approx u_n \sum_{m=0}^{N-1}v(m)W^{-nm}_{N} \\ & \approx u_nv_n \\ ...\\ DFT(u_nv_n) &= \sum_{n=0}^{N-1}u_nv_nW^{kn}_{N} \approx \sum_{m=0}^{N-1} U(k-m)V(m) \end{aligned}
    与上一节推导类似 ,这里 k=mm+N1U(km)WNn(km)\sum_{k=m}^{m+N-1} U(k-m)W^{-n(k-m)}_{N}也是抽取 k=0,1,2,…N-1-m 的频普(N-m个),补充m个零幅值频率点再做傅里叶逆变换,所以也就是约等于unu_n(一些高频频率丢失),当卷积点数远小于频率数量时,这个约等于是成立的(V(k)也要在后面填充N-M个零)。

    结论:对两组离散数据相乘,相当于其两组数据频域作卷积,也就是起到频域数据滤波的效果

    DFT与窗函数
    由上一个公式可知道,离散数字的乘积与傅里变换后卷积的关系,
    unvnConvolve(U(k),V(k))u_nv_n \approx Convolve(U(k),V(k))
    窗函数v与原始数据u相乘后为x,x的频域数据相当于u的频域数据做了一个卷积后的结果


    用python 做一下实验

    import numpy as np
    import math
    %pylab inline
    
    data = np.random.randn(64)  #尺寸应该为2的次方
    kernel = np.array([0.2,0.5,0.2])
    convres = np.convolve(data,kernel,'same') #卷积后的结果
    
    dataft = np.fft.fft(data) #原始数据傅里叶结果
    kernelft = np.fft.fft(np.pad(kernel,(0,data.size-kernel.size),'constant')) #卷积核傅里叶结果 (填充0 使其与原始数据长度一致)
    convresft = np.fft.fft(convres) #卷积后数据傅里叶结果
    approx_times = kernelft*dataft #原始数据傅里叶变换与卷积核傅里叶变换相乘结果
    
    div = convresft/(approx_times+1e-10)
    div:np.ndarray
    #可以看到 kernelft*dataft 与 convresft是比较接近的,(幅值很相近 可能相位相差大一点)
    print('可以看到 kernelft*dataft 与 convresft是比较接近的,(幅值很相近 可能相位相差大一点)')
    print('幅值比 = ',np.abs(div)) 
    print('相位比 = ',np.angle(div))
    
    可以看到 kernelft*dataft 与 convresft是比较接近的,(幅值很相近 可能相位相差大一点)
    幅值比 =  [ 1.08514358  1.10805967  0.95296099  1.03198588  0.87026809  0.99448745
      1.00526948  0.95799491  1.05421725  1.29726256  0.80609292  0.88131322
      0.91631737  1.03311114  0.82737111  0.86807867  1.42545673  0.56157421
      0.91429024  1.11167318  0.91778763  1.11124952  1.14813349  0.95219242
      1.05733615  0.97868321  1.23478938  1.42048091  1.14337804  2.09690390
      2.00851084 12.23768621  1.64446639 12.23768621  2.00851084  2.09690390
      1.14337804  1.42048091  1.23478938  0.97868321  1.05733615  0.95219242
      1.14813349  1.11124952  0.91778763  1.11167318  0.91429024  0.56157421
      1.42545673  0.86807867  0.82737111  1.03311114  0.91631737  0.88131322
      0.80609292  1.29726256  1.05421725  0.95799491  1.00526948  0.99448745
      0.87026809  1.03198588  0.95296099  1.10805967]
    相位比 =  [ 0.00000000  0.20103405  0.28134838  0.21410434  0.41299861  0.41863674
      0.73898843  0.64316977  0.74492692  0.78055698  1.08688623  0.98047538
      1.14141382  1.17736177  1.37149334  1.44358457  1.38108547  1.87904898
      1.75100384  1.82650733  1.70435237  1.78623087  2.11383998  1.87868604
      2.14499912  2.06628397  2.81174815  3.10020892  2.29021613 -2.58993841
      3.12919688 -1.20352958  3.14159265  1.20352958 -3.12919688  2.58993841
     -2.29021613 -3.10020892 -2.81174815 -2.06628397 -2.14499912 -1.87868604
     -2.11383998 -1.78623087 -1.70435237 -1.82650733 -1.75100384 -1.87904898
     -1.38108547 -1.44358457 -1.37149334 -1.17736177 -1.14141382 -0.98047538
     -1.08688623 -0.78055698 -0.74492692 -0.64316977 -0.73898843 -0.41863674
     -0.41299861 -0.21410434 -0.28134838 -0.20103405]
    
    windows = np.ndarray(data.size)
    for i in range(windows.size):
        windows[i] = 0.5*(1-math.cos(2*math.pi*i/windows.size)) #汉宁窗
        
    dataft = np.fft.fft(data)/data.size #原始数据傅里叶结果 除以尺寸是因为FFT结果是频域幅的N倍
    windowsft = np.fft.fft(windows)/windows.size #窗数据傅里叶结果
    #windowsft = np.concatenate([windowsft[0:2],windowsft[-2:-1]])# windowsft[0:2] #由于汉宁窗只有前两个频率是有效的
    
    windata = windows * data #经过窗函数后的数据
    windataft = np.fft.fft(windata)/windata.size #开窗后 傅里叶变换数据
    
    #conv1res = np.convolve(dataft,windowsft,'same')
    #conv1res = np.convolve(dataft,windowsft,'full')[windowsft.size-1:windowsft.size-1+64] #原数据傅里叶变换与窗函数傅里叶变换卷积
    
    conv1res = np.convolve(dataft,windowsft,'full')[0:dataft.size]
    
    div1 = windataft/conv1res
    np.set_printoptions(suppress=True)
    np.set_printoptions(precision=5)
    print('可以看到 windataft 与 conv1res是比较接近的,(幅值很大部分频率位置相近 可能相位相差大一点)')
    print('幅值比 = ',np.abs(div1))
    print('相位比 = ',np.angle(div1))
    
    
    可以看到 windataft 与 conv1res是比较接近的,(幅值很大部分频率位置相近 可能相位相差大一点)
    幅值比 =  [0.58146 1.51554 1.28522 0.97078 1.56031 1.37292 0.97845 1.65274 1.1212
     1.66299 0.64529 0.81242 1.39261 1.29182 0.31025 1.21369 1.08077 0.97666
     0.57983 1.86432 1.35158 0.82371 4.23688 1.03864 0.78796 0.99415 2.48601
     0.7744  1.50276 1.1419  1.02537 0.96125 0.43243 1.09725 1.29529 0.87585
     1.30709 1.07283 1.27008 0.36948 1.30446 1.0443  0.71414 1.50539 1.12606
     0.97264 1.1533  1.08976 2.0767  0.87196 0.65309 1.43666 0.71085 0.73239
     1.17037 1.56127 1.22908 1.19577 1.43997 1.2673  1.01911 1.31521 1.40685
     1.     ]
    相位比 =  [ 0.      -1.61908 -0.08994 -0.38359 -0.1287  -0.04803 -0.17619 -0.10482
     -0.04996  0.45878 -0.67394  0.36607 -0.01151 -0.14093 -1.83447 -0.02431
      0.20923  0.17784  0.15037 -0.23009 -0.23083 -0.03418 -2.95075  0.26942
     -0.09524  0.59162 -0.55638  0.39722 -0.93279  0.08441  0.56759 -0.13494
      1.68964  0.15164  0.75175 -0.26193  0.60153 -0.07737 -0.03958  0.63195
      1.05636 -0.13213  0.48739 -0.25624 -0.16226  0.20864  1.42779  0.40813
      0.18961 -0.03562 -0.66101 -0.47631  0.97983 -0.46818  0.10568  0.13711
      0.09992 -0.47046 -0.14611  0.00341 -0.275   -0.31533  0.0306   0.     ]
    
    wewe = np.fft.ifft(conv1res*data.size)
    plt.subplot(311)
    plt.plot(data)    #原始数据波形
    plt.subplot(312)  
    plt.plot(windata)  #乘以窗数据后的波形
    plt.subplot(313)
    plt.plot(np.real(wewe))   #傅里叶结果卷积后 再傳里叶逆变换的波形
    

    在这里插入图片描述

    展开全文
  • cvxpy工具包求解稀疏约束优化问题 目录 快速傅里叶变换(FFT) 窗函数傅里叶变换的影响 两种信号滤波方式 短时傅里叶变换(STFT) 1. 快速傅里叶变换 FFT的实现主要有两种方法,分别通过np.fft或scipy.fftpack中的fft...

    摘要:一直以来都是用MATLAB做信号处理,得到预处理的特征后再用Python进一步应用神经网络之类的方法。这里将MATLAB中的FFT、STFT、加窗以及带通滤波通过Python接口实现,防止以后MATLAB用不了了,一定程度上也提高了效率,不用两个软件换来换去。

    系列目录

    1. Python信号处理:快速傅里叶变换(FFT),短时傅里叶变换(STFT),窗函数,以及滤波
    2. Python信号处理:自相关函数(对标MATLAB中的autocorr)
    3. Python信号处理:波束形成及目标方位估计,CBF、MVDR
    4. Python信号处理:cvxpy工具包求解稀疏约束优化问题

    目录

    1. 快速傅里叶变换(FFT)
    2. 窗函数对傅里叶变换的影响
    3. 两种信号滤波方式
    4. 短时傅里叶变换(STFT)

    1. 快速傅里叶变换

    FFT的实现主要有两种方法,分别通过np.fftscipy.fftpack中的fft实现,实现方法基本是一样的,以scipy.fftpack.fft为例。

    首先生成一个带噪声的信号,频率为10、15、20 Hz,采样率fsf_s=100 Hz,时长2 s。时域波形如下左图所示。使用scipy.fftpack.fft获取频域信号,同时使用scipy.fftpack.fftfreq获取,信号频谱如下右图所示。

    X = fftpack.fft(x)           # x为时域信号
    freqs = fftpack.fftfreq(N) * fs  # N为信号点数,fs为采样率
    

    在这里插入图片描述

    2. 窗函数对傅里叶变换的影响

    以凯撒窗np.kaiser(点数,系数)为例,凯撒窗能抑制旁瓣,但会加宽主瓣。使用时,在fft中对信号加窗即可。下图给出不同系数的凯撒窗对主瓣和旁瓣的影响。可以看到,加窗后旁瓣受到抑制,但主瓣变宽。通过调整系数,可以改变凯撒窗的窗口。

    X = fftpack.fft(np.kaiser(N, 6) * x)
    

    在这里插入图片描述

    3. 两种信号滤波方式

    A. 通过FFT和IFFT进行滤波

    先对信号进行FFT,然后保留想要的频率,将不要的频率都设为0,再将频谱进行IFFT得到滤波后的时域波形。

    def filter1(signal, freq, freqs):
        """利用FFT进行滤波。
        # Arguments
        #   signal: 带噪声的时域信号
        #   freq: 想要的频率
        #   freqs: 信号频率
        # Return
        #   y: 滤波后的时域信号
        """
        X = fftpack.fft(signal)
        index_freq = [np.where(i == freqs)[0][0] for i in freq]
        index_zeros = np.setdiff1d(np.arange(0, len(freqs)), index_freq)
        X[index_zeros] = 0
        y = fftpack.ifft(X)
        return y
    
    y = filter1(x, [10, 15, 20], freqs)
    

    滤波后的时域信号及其频谱如下图所示。
    在这里插入图片描述

    B. 通过butter和filtfilt进行滤波

    使用scipy.signal.butter通过截至频率获取滤波器系数,利用scipy.signal.filtfilt进行滤波。假设频率为10、15、20 Hz,采样率fsf_s=100 Hz,时长2 s的时域信号,截止频率为10和20 Hz。

    wl = 2 * (10 - 1) / fs
    wh = 2 * (20 + 1) / fs
    b, a = signal.butter(8, [wl ,wh], 'bandpass')  #8表示滤波器的阶数
    filtedData = signal.filtfilt(b, a, x)
    

    滤波后的时域信号及其频谱如下图所示。
    在这里插入图片描述

    4. 短时傅里叶变换

    scipy.signal.stftscipy.signal.spectrogram两个类似的方法,区别在于后者返回的是能量,前者返回的是频谱,差一个绝对值的平方。

    两者接口参数如下:

    • scipy.signal.stft(x,fs=1.0,window='hann',nperseg=256,noverlap=None,nfft=None,detrend=False,return_onesided=True,boundary ='zeros',padded=True,axis=-1 )

    • scipy.signal.spectrogram(x, fs=1.0, window=('tukey', 0.25), nperseg=None, noverlap=None, nfft=None, detrend=’constant’, return_onesided=True, scaling=’density’, axis=-1, mode=’psd’)

    常用参数从左往右依次为:信号、采样率、窗、每段长度、重叠长度、FFT点数、消除趋势,两者基本都是一样的。

    展开全文
  • 2、快速傅里叶变换 3、信号窗函数 4、卷积 导入库 import matplotlib.pyplot as plt import scipy.signal as sgn import numpy as np 1、过滤:以某种方式修改输入信号 1)快速线性两次应用滤波函数 filtfilt() ...
  • 文章目录傅里叶变换及傅里叶逆变换定义窗函数/矩形脉冲信号的傅里叶变换基于MATLAB的快速傅里叶变换FFT 傅里叶变换及傅里叶逆变换定义 能从时域的非周期连续信号转化到频域非周期连续信号。 窗函数/矩形脉冲...
  • 为了改善窗口傅里叶变换滤波算法中的阈值选取方法,提出了自适应阈值方法。通过模拟散斑干涉相位图验证可知,所提出的自适应阈值方法可以对噪声进行有效的滤除,同时对信号信息进行很好的保留。
  • 对分析方法加上窗函数,只分析在窗内的信号,大致上可以说是实现了定位,但若想要实现精确的定位,必须能随时改变时间窗的大小,短时傅里叶变换是无法做到的,因此就产生了小波理论、多分辨分析,其可以实现对信号在...
  • \qquad设g(t)=g(−t)g(t)=g(-t)g(t)=g(−t)是一个实的对称窗函数,对它平移uuu各单位,以及调制(7-1小波分析之傅里叶分析有介绍平移与调制等性质)得到gu,ξ(t)=e(iξt)g(t−u)g_{u,\xi}(t)=e^{(i\xi t)}g(t-u)gu,...
  • MATLAB实现短时傅里叶变换

    千次阅读 2020-10-19 15:54:37
    一、短时傅里叶变换的定义 离散傅里叶变换使用的是一种全局变换,因为它表示一段时间内平均的频率特性,无法表述信号的时域局域性质,为了能够分析处理非平稳...1.短时傅里叶变换tfrstft函数: [tfr,t,f]=tfrstft(x,

空空如也

空空如也

1 2 3 4 5 ... 14
收藏数 262
精华内容 104
关键字:

傅里叶变换窗函数