zoukankan      html  css  js  c++  java
  • SciPy fftpack(傅里叶变换)

    章节


    SciPy提供了fftpack模块,包含了傅里叶变换的算法实现。

    傅立叶原理表明:任何连续测量的时序或信号,都可以表示为不同频率的正弦波信号的无限叠加。傅里叶变换把信号从时域变换到频域,以便对信号进行处理。傅里叶变换在信号与噪声处理、图像处理、音频信号处理等领域得到了广泛应用。

    如需进一步了解傅里叶变换原理,可以参考相关资料。

    快速傅里叶变换

    计算机只能处理离散信号,使用离散傅里叶变换(DFT) 是计算机分析信号的基本方法。但是离散傅里叶变换的缺点是:计算量大,时间复杂度太高,当采样点数太高的时候,计算缓慢,由此出现了DFT的快速实现,即快速傅里叶变换FFT。

    快速傅里叶变换(FFT)是计算量更小的离散傅里叶变换的一种实现方法,其逆变换被称为快速傅里叶逆变换(IFFT)。

    示例

    print(fft(np.array([4., 3., 5., 10., 5., 3.])))

    先对数据进行fft变换,然后再ifft逆变换。

    import numpy as np
    #从fftpack中导入fft(快速傅里叶变化)和ifft(快速傅里叶逆变换)函数
    from scipy.fftpack import fft,ifft
    
    #创建一个随机值数组
    x = np.array([1.0, 2.0, 1.0, -1.0, 1.5])
    
    #对数组数据进行傅里叶变换
    y = fft(x)
    print('fft: ')
    print(y)
    print('
    ')
    
    #快速傅里叶逆变换
    yinv = ifft(y)
    print('ifft: ')
    print(yinv)
    print('
    ')
    

    输出

    
    fft:
    [ 4.5       +0.j          2.08155948-1.65109876j -1.83155948+1.60822041j
     -1.83155948-1.60822041j  2.08155948+1.65109876j]
    
    
    ifft:
    [ 1. +0.j  2. +0.j  1. +0.j -1. +0.j  1.5+0.j]
    
    

    可以看到fft,ifft返回的都是复数。ifft返回的结果中,复数的虚部都是0,实部与原始数据x一致。

    这些点的频率无法计算,因为没有设置这N个点的时间长度。如不理解,不必深究,后面会介绍。

    理解fft变换结果

    我们知道,傅里叶变换把时域信号变为频域信号。在离散傅里叶变换中,频域信号由一系列不同频率的谐波(频率成倍数)组成。fft返回值是一个复数数组,每个复数表示一个正弦波。通常一个波形由振幅,相位,频率三个变量确定,可以从fft的返回值里,获取这些信息。

    假设a是时域中的周期信号,采样频率为Fs,采样点数为N。如果A[N] = fft(a[N]),返回值A[N]是一个复数数组,其中:

    • A[0]表示频率为0hz的信号,即直流分量。
    • A[1:N/2]包含正频率项,A[N/2:]包含负频率项。正频率项就是转化后的频域信号,通常我们只需要正频率项,即前面的n/2项,负频率项是计算的中间结果(正频率项的镜像值)。
    • 每一项的频率计算:假设A[i]为数组中的元素,表示一个波形,该波形的频率 = i * Fs / N
    • A[i] = real + j * imag,是一个复数,相位就是复数的辐角,相位 = arg(real/imag)
    • 类似的,振幅就是复数的模,振幅 = sqrt(real2+imag2)。但是fft的返回值的模是放大值,直流分量的振幅放大了N倍,弦波分量的振幅放大了N/2倍。

    频率分辨率
    频率分辨率是离散傅里叶变换(DFT)频域相邻刻度之间的实际频率之差。采样时,数据采样了T秒(T = 采样点数N / 采样频率Fs),信号的成分中周期最大也就是T秒,最低频率即“基频”就等于1 / T,也就是Fs / N,这就是频率分辨率。基频 = Fs / N,各个谐波的频率就是 i * Fs / N,这个公式用于计算各个波形的频率。

    示例

    import numpy as np
    from scipy.fftpack import fft
    
    # 采样点数
    N = 4000
    
    # 采样频率 (根据采样定理,采样频率必须大于信号最高频率的2倍,信号才不会失真)
    Fs = 8000
    x = np.linspace(0.0, N/Fs, N)
    
    # 时域信号,包含:直流分量振幅1.0,正弦波分量频率100hz/振幅2.0, 正弦波分量频率150Hz/振幅0.5/相位np.pi
    y = 1.0 + 2.0 * np.sin(100.0 * 2.0*np.pi*x) + 0.5*np.sin(150.0 * 2.0*np.pi*x + np.pi)
    
    # 进行fft变换
    yf = fft(y)
    
    # 获取振幅,取复数的绝对值,即复数的模
    abs_yf = np.abs(yf)
    
    # 获取相位,取复数的角度
    angle_y=np.angle(yf)
    
    # 直流信号
    print('
    直流信号')
    print('振幅:', abs_yf[0]/N) # 直流分量的振幅放大了N倍
    
    # 100hz信号
    index_100hz = 100 * N // Fs # 波形的频率 = i * Fs / N,倒推计算索引:i = 波形频率 * N / Fs
    print('
    100hz波形')
    print('振幅:', abs_yf[index_100hz] * 2.0/N) # 弦波分量的振幅放大了N/2倍
    print('相位:', angle_y[index_100hz])
    
    # 150hz信号
    index_150hz = 150 * N // Fs # 波形的频率 = i * Fs / N,倒推计算索引:i = 波形频率 * N / Fs
    print('
    150hz波形')
    print('振幅:', abs_yf[index_150hz] * 2.0/N) # 弦波分量的振幅放大了N/2倍
    print('相位:', angle_y[index_150hz])
    print('100hz与150hz相位差:',  angle_y[index_150hz] - angle_y[index_100hz])
    print('
    ')
    
    
    

    输出

    
    直流信号
    振幅: 1.0
    
    100hz波形
    振幅: 1.9989359813189005
    相位: -1.5315264186250062
    
    150hz波形
    振幅: 0.5008489983048182
    相位: 1.6297011890497097
    100hz与150hz相位差: 3.161227607674716
    
    

    可以看到,正弦波的相位不一定从0开始,但波形之间的相位差确实s约等于一个pi(值跟采样频率与采样点数有关系)。

    离散余弦变换(DCT)

    由于许多要处理的信号都是实信号,在使用FFT时,对于实信号,傅立叶变换的共轭对称性导致在频域中有一半的数据冗余。

    离散余弦变换(DCT)是对实信号定义的一种变换,变换后在频域中得到的也是一个实信号,相比离散傅里叶变换DFT而言, DCT可以减少一半以上的计算。DCT还有一个很重要的性质(能量集中特性):大多书自然信号(声音、图像)的能量都集中在离散余弦变换后的低频部分,因而DCT在(声音、图像)数据压缩中得到了广泛的使用。由于DCT是从DFT推导出来的另一种变换,因此许多DFT的属性在DCT中仍然是保留下来的。

    SciPy.fftpack中,提供了离散余弦变换(DCT)与离散余弦逆变换(IDCT)的实现。

    示例

    import numpy as np
    from scipy.fftpack import dct,idct
    y = dct(np.array([4., 3., 5., 10., 5., 3.]))
    print(y)
    

    输出

    [ 60.          -3.48476592 -13.85640646  11.3137085    6.
      -6.31319305]
    

    离散余弦逆变换(idct),是离散余弦变换(DCT)的反变换。

    示例

    import numpy as np
    from scipy.fftpack import dct,idct
    y = idct(np.array([4., 3., 5., 10., 5., 3.]))
    print(y)
    

    输出

    [ 39.15085889 -20.14213562  -6.45392043   7.13341236   8.14213562
      -3.83035081]
    
  • 相关阅读:

    《银行》
    动画
    定位
    浮动
    《css的总结》
    Java的Debug调试
    Eclipse "Unable to install breakpoint due to missing line number attributes..."
    解决 Tomcat Server in Eclipse unable to start within 45 seconds 不能启动的问题
    Java下载文件(流的形式)
  • 原文地址:https://www.cnblogs.com/jinbuqi/p/11803923.html
Copyright © 2011-2022 走看看