zoukankan      html  css  js  c++  java
  • 快速傅里叶变换模块(fft)

    什么是傅里叶变换?

    法国科学家傅里叶提出,任何一条周期曲线,无论多么跳跃或不规则,都能表示成一组光滑正弦曲线叠加之和。

    傅里叶变换的目的是可将时域(即时间域)上的信号转变为频域(即频率域)上的信号,随着域的不同,对同一个事物的了解角度也就随之改变,因此在时域中某些不好处理的地方,在频域就可以较为简单的处理。这就可以大量减少处理信号存储量。

    例如:弹钢琴

    假设有一时间域函数:y = f(x),根据傅里叶的理论它可以被分解为一系列正弦函数的叠加,他们的振幅A,频率ω或初相位φ不同:

    所以傅里叶变换可以把一个比较复杂的函数转换为多个简单函数的叠加,看问题的角度也从时间域转到了频率域,有些的问题处理起来就会比较简单。

    傅里叶变换相关函数

    导入快速傅里叶变换所需模块

    import numpy.fft as nf

    通过采样数与采样周期求得傅里叶变换分解所得曲线的频率序列

    freqs = np.fft.fftfreq(采样数量, 采样周期)

    通过原函数值的序列j经过快速傅里叶变换得到一个复数数组,复数的模代表的是振幅,复数的辐角代表初相位

    np.fft.fft(原函数值序列) -> 目标函数值序列(复数)

    通过一个复数数组(复数的模代表的是振幅,复数的辐角代表初相位)经过逆向傅里叶变换得到合成的函数值数组

    np.fft.ifft(目标函数值序列(复数))->原函数值序列

    案例:针对合成波做快速傅里叶变换,得到一组复数序列;再针对该复数序列做逆向傅里叶变换得到新的合成波并绘制。

    ffts = nf.fft(sigs6)
    sigs7 = nf.ifft(ffts).real
    mp.plot(times, sigs7, label=r'$omega$='+str(round(1 / (2 * np.pi),3)), alpha=0.5, linewidth=6)

    案例:针对合成波做快速傅里叶变换,得到分解波数组的频率、振幅、初相位数组,并绘制频域图像。

    # 得到分解波的频率序列
    freqs = nf.fftfreq(times.size, times[1] - times[0])
    # 复数的模为信号的振幅(能量大小)
    ffts = nf.fft(sigs6)
    pows = np.abs(ffts)
    
    mp.subplot(122)
    mp.title('Frequency Domain', fontsize=16)
    mp.xlabel('Frequency', fontsize=12)
    mp.ylabel('Power', fontsize=12)
    mp.tick_params(labelsize=10)
    mp.grid(linestyle=':')
    mp.plot(freqs[freqs >= 0], pows[freqs >= 0], c='orangered', label='Frequency Spectrum')
    mp.legend()
    mp.tight_layout()
    mp.show()

     

    """
    傅里叶变换拆方波
    """
    import numpy as np
    import matplotlib.pyplot as mp
    import numpy.fft as nf
    
    x = np.linspace(-2*np.pi, 2*np.pi, 1000)
    y = np.zeros(x.size)
    for i in range(1, 1000):
        y += 4*np.pi/(2*i-1) * np.sin((2*i-1)*x)
    
    mp.figure('FFT',facecolor='lightgray')
    mp.subplot(121)
    mp.title('Time Domain',fontsize=18)
    mp.grid(linestyle=":")
    mp.plot(x, y, label=r'$y$')
    #针对方波y,做fft
    comp_ary = nf.fft(y)
    y2 = nf.ifft(comp_ary)
    mp.plot(x,y2,color = 'orangered',linewidth=5,alpha=0.5,label='$y$')
    #绘制频域图像(频率/能量)
    mp.subplot(122)
    freqs = nf.fftfreq(y.size,x[1]-x[0])
    pows = np.abs(comp_ary)#负数的模
    mp.title('Frequency Domain',fontsize=18)
    mp.grid(linestyle=":")
    mp.plot(freqs[freqs>0],pows[freqs>0],color= 'orangered',label='frequency')
    
    mp.legend()
    mp.show()

    基于傅里叶变换的频域滤波

    含噪信号是高能信号与低能噪声叠加的信号,可以通过傅里叶变换的频域滤波实现降噪。

    通过FFT使含噪信号转换为含噪频谱,去除低能噪声,留下高能频谱后再通过IFFT留下高能信号。

    案例:基于傅里叶变换的频域滤波为音频文件去除噪声。

    读取音频文件,获取音频文件基本信息:采样个数,采样周期,与每个采样的声音信号值。绘制音频时域的:时间/位移图像。

    import numpy as np
    import numpy.fft as nf
    import scipy.io.wavfile as wf
    import matplotlib.pyplot as mp
    
    sample_rate, noised_sigs = wf.read('../data/noised.wav')
    noised_sigs = noised_sigs / 2 ** 15
    times = np.arange(len(noised_sigs)) / sample_rate
    mp.figure('Filter', facecolor='lightgray')
    mp.subplot(221)
    mp.title('Time Domain', fontsize=16)
    mp.ylabel('Signal', fontsize=12)
    mp.tick_params(labelsize=10)
    mp.grid(linestyle=':')
    mp.plot(times[:178], noised_sigs[:178],c='orangered', label='Noised')
    mp.legend()
    mp.show()

    基于傅里叶变换,获取音频频域信息,绘制音频频域的:频率/能量图像。

    freqs = nf.fftfreq(times.size, 1 / sample_rate)
    noised_ffts = nf.fft(noised_sigs)
    noised_pows = np.abs(noised_ffts)
    mp.subplot(222)
    mp.title('Frequency Domain', fontsize=16)
    mp.ylabel('Power', fontsize=12)
    mp.tick_params(labelsize=10)
    mp.grid(linestyle=':')
    mp.semilogy(freqs[freqs >= 0],noised_pows[freqs >= 0], c='limegreen',label='Noised')
    mp.legend()

    将低能噪声去除后绘制音频频域的:频率/能量图像。

    fund_freq = freqs[noised_pows.argmax()]
    noised_indices = np.where(freqs != fund_freq)
    filter_ffts = noised_ffts.copy()
    filter_ffts[noised_indices] = 0
    filter_pows = np.abs(filter_ffts)
    
    mp.subplot(224)
    mp.xlabel('Frequency', fontsize=12)
    mp.ylabel('Power', fontsize=12)
    mp.tick_params(labelsize=10)
    mp.grid(linestyle=':')
    mp.plot(freqs[freqs >= 0], filter_pows[freqs >= 0],c='dodgerblue', label='Filter')
    mp.legend() 

    基于逆向傅里叶变换,生成新的音频信号,绘制音频时域的:时间/位移图像。

    filter_sigs = nf.ifft(filter_ffts).real
    mp.subplot(223)
    mp.xlabel('Time', fontsize=12)
    mp.ylabel('Signal', fontsize=12)
    mp.tick_params(labelsize=10)
    mp.grid(linestyle=':')
    mp.plot(times[:178], filter_sigs[:178],c='hotpink', label='Filter')
    mp.legend()

    重新生成音频文件。

    wf.write('../../data/filter.wav',sample_rate,(filter_sigs * 2 ** 15).astype(np.int16))
    """
    频域滤波
    """
    import numpy as np
    import numpy.fft as nf
    import matplotlib.pyplot as mp
    import scipy.io.wavfile as wf
    
    # 采样率,采样位移
    sample_rate, sigs = wf.read('noised.wav')
    print(sample_rate)  # 44100
    print(sigs.shape)  # (220500,)
    sigs = sigs/(2**15)
    times = np.arange(sigs.size) / sample_rate
    print(times)
    """
    [0.00000000e+00 2.26757370e-05 4.53514739e-05 ... 4.99993197e+00
     4.99995465e+00 4.99997732e+00]
    """
    # 绘图
    mp.figure('Filter', facecolor='lightgray')
    mp.subplot(221)
    mp.title('Time Domain', fontsize=16)
    mp.ylabel('Signal', fontsize=14)
    mp.grid(linestyle=':')
    mp.plot(times[:178], sigs[:178], color='dodgerblue', label='Noised Signal')
    mp.legend()
    
    # fft变换后绘制频域图(频率/能量)
    comp_ary = nf.fft(sigs)
    pows = np.abs(comp_ary)
    freqs = nf.fftfreq(sigs.size, times[1] - times[0])
    mp.subplot(222)
    mp.title('Frequency Domain', fontsize=16)
    mp.ylabel('pow', fontsize=14)
    mp.grid(linestyle=":")
    mp.semilogy(freqs[freqs > 0], pows[freqs > 0], color='orangered', label='Freq')
    
    mp.legend()
    
    # 在频谱中去除噪声
    maxpow_freq = freqs[pows.argmax()]
    noised_inds = np.where(freqs != maxpow_freq)
    # 噪声的复数全抹为0
    comp_ary[noised_inds] = 0
    # 绘制降噪之后的频谱图像
    pows = np.abs(comp_ary)
    mp.subplot(224)
    mp.ylabel('pow', fontsize=14)
    mp.grid(linestyle=':')
    mp.plot(freqs[freqs > 0], pows[freqs > 0],
            color='orangered', label='Freq')
    mp.legend()
    
    
    # 逆向傅里叶变换,输出时域函数图像
    filter_sigs = nf.ifft(comp_ary)
    mp.subplot(223)
    mp.ylabel('Signal', fontsize=14)
    mp.grid(linestyle=":")
    mp.plot(times[:178], filter_sigs[:178], color='dodgerblue', label='Filter Signal')
    mp.legend()
    # 保存文件
    wf.write('filter.wav', sample_rate, (filter_sigs * 2 ** 15).astype(np.int16))
    
    mp.tight_layout()
    mp.show()

     

  • 相关阅读:
    AutoCAD LISP矩形窗格绘制
    AutoCAD VBA多重延伸
    2011年3月24日星期四
    AutoCAD VBA对齐对象
    AutoCAD VBA根据对象缩放
    AutoCAD LISP绘制多个等半径圆相切
    AutoCAD LISP利用一顶点和三边长绘制三角形
    AutoCAD VBA对象的组合和拆散
    AutoCAD LISP修改已存在圆半径
    AutoCAD VBA基于对象的分层
  • 原文地址:https://www.cnblogs.com/maplethefox/p/11490143.html
Copyright © 2011-2022 走看看