一、基于傅里葉定理，用一組正弦函數合成方波

　　

'''
    三角函數通用函數
    傅里葉定理：任何一個周期性曲線，無論多么跳躍或者不規則，都可以被解析成一組光滑的正弦函數的疊加
            ---應用：合成方波(即不規則的方波由一組光滑的正弦函數疊加合成的)
                    如：y = 4π/(2*n-1) * sin((2*n-1)*x)
'''
import numpy as np
import matplotlib.pyplot as mp

x = np.linspace(-2 * np.pi, 2 * np.pi, 1000)
y1 = 4 * np.pi * np.sin(x)
y2 = 4 * np.pi / 3 * np.sin(3 * x)

n = 6
y = np.zeros(1000)
for i in range(1, n + 1):
    y += 4 * np.pi / (2 * i - 1) * np.sin((2 * i - 1) * x)

mp.plot(x, y1, label='y1', alpha=0.3)
mp.plot(x, y2, label='y2', alpha=0.3)
mp.plot(x, y, label='y')
mp.legend()
mp.show()

　　

　　

二、基于快速傅里葉變換實現方波的拆解

　　

'''
    傅里葉變換：即是將一條周期性忑基于傅里葉定理進行拆解，得到一組光滑的正弦曲線的變換過程
    傅里葉變化的目的：是可將時間域的信號轉換為頻域（頻率域）信號，隨著域的不同，對同一個事物的了解角度也就隨之改變，
                    因此在時域中某些不好處理的地方，在頻域中可以較為簡單的處理，如：數據存儲、數據降噪等
    傅里葉變換主要運用在音頻文件領域
    在numpy中有一個專門模塊---fft模塊，專門處理傅里葉變換
        ---復數數組 = fft.fft(原函數y的數組)  --->即為快速傅里葉變換
            --復數數組中的每個復數可以描述一條正弦曲線，復數的模代表振幅A，復數的輔角代表初相位角φ
        ---freqs = fft.fftfreq(采樣數量，采樣周期) -->通過采樣數和采樣周期求得曲線的頻率序列，即ω

        ---逆向傅里葉變換
            -- y2 = fft.ifft(復數數組)

    案例：拆解方波，繪制圖像
'''

import numpy as np
import matplotlib.pyplot as mp
import numpy.fft as fft

x = np.linspace(-2 * np.pi, 2 * np.pi, 1000)

# 合成方波
n = 1000
y = np.zeros(1000)
for i in range(1, n + 1):
    y += 4 * np.pi / (2 * i - 1) * np.sin((2 * i - 1) * x)

# 拆解方波
freqs = fft.fftfreq(y.size, y[1] - y[0])
complex_ary = fft.fft(y)
print(complex_ary.size)
# 逆向傅里葉變換
y2 = fft.ifft(complex_ary)

# 繪制時域圖
mp.figure('FFT')
mp.subplot(121)
mp.title('Time Domain')
mp.plot(x, y, label='y')
mp.plot(x, y2, linewidth=5, label='y2', alpha=0.4)
mp.grid(linestyle=":")
mp.legend()

# 繪制頻域圖
mp.subplot(122)
mp.title('Frequency Domain')
mp.grid(linestyle=":")
power = np.abs(complex_ary)
mp.plot(freqs[freqs > 0], power[freqs > 0])

mp.tight_layout()
mp.show()

　　

三、基于傅里葉變換的頻域濾波

　　

'''
    基于傅里葉變換的頻域濾波
        ----過濾音頻中噪聲：含噪信號由高能信號與低能噪聲疊加而成，可以通過FFT的頻域濾波實現降噪。通過FFT使含噪信號轉為含噪頻譜，
            留下高能頻譜，過濾掉低能噪聲，再經過IFFT將高能頻譜生成高能信號。
    步驟：
        1.讀取音頻文件，獲取音頻文件的基本信息：采樣率、采樣周期、每個采樣點的聲音位移值。繪制音頻的時域（時間/位移）圖像
        2.基于傅里葉變換，獲取音頻頻域信息，繪制音頻的頻域（頻率/能量）圖像
        3.將低能噪聲去除，繪制音頻頻域（頻域/能量）圖像
        4.基于逆向傅里葉變換（ifft)，生成新的音頻信號，繪制圖像
        5.生成音頻文件

'''
import scipy.misc as sc
import matplotlib.pyplot as mp
import numpy as np
import scipy.io.wavfile as wf  # 讀取音頻文件模塊
import numpy.fft as nf  # 傅里葉變換模塊
import warnings

warnings.filterwarnings("ignore")

# 第一步驟
# 讀取音頻文件，獲取采樣率(每秒采樣多少個點）和采樣點位移（一共采樣多少個點）
sample_rate, noise_sigs = wf.read('./da_data/noised.wav')
# print('sample_rate:', sample_rate)
# print('noise_sigs:', noise_sigs.shape)
# 計算采樣時間
times = np.arange(len(noise_sigs)) / sample_rate
print(times.shape)

mp.figure('Filter', facecolor='lightgray')
mp.subplot(221)
mp.title('Time Domain')
mp.xlabel('Time')
mp.ylabel('Signal')
mp.grid(linestyle=':')
# 因數據太大，只取前178個
mp.plot(times[:178], noise_sigs[:178], c='dodgerblue', label='Signals')
mp.legend()

# 第二步驟
freqs = nf.fftfreq(times.size, 1 / sample_rate)
complex_ary = nf.fft(noise_sigs)
# 求模，代表能量大小
powers = np.abs(complex_ary)
mp.subplot(222)
mp.title('Frequent Domain')
mp.ylabel('power')
mp.grid(linestyle=":")
mp.semilogy(freqs[freqs > 0], powers[freqs > 0], c='orangered', label='Noised')
mp.legend()

# 第三步驟
# 找到頻率最大值
fund_freq = freqs[powers.argmax()]
# 找出所有噪聲下標,where方法返回的是下標數組
nosie_indexes = np.where(freqs != fund_freq)
# 拷貝一份，避免修改原始數據
complex_filter = complex_ary.copy()
# 將噪聲賦值0---濾波
complex_filter[nosie_indexes] = 0
# 求濾波的模大小,代表能量大小
power_filter = np.abs(complex_filter)
# 畫圖
mp.subplot(224)
mp.ylabel('power')
mp.grid(linestyle=":")
mp.plot(freqs[freqs > 0], power_filter[freqs > 0], c='orangered', label='Filter')
mp.legend()

# 第四步驟
# 根據反向傅里葉，得到過濾后的采樣點位移
filter_sigs = nf.ifft(complex_filter).real
# 畫圖
mp.subplot(223)
mp.ylabel('Signal')
mp.grid(linestyle=':')
# 因數據太大，只取前178個
mp.plot(times[:178], filter_sigs[:178], c='dodgerblue', label='Signals')
mp.legend()

# 第五步驟
# 參數需要，保存路徑、采樣率、采樣位移
wf.write('./da_data/out.wav', sample_rate, filter_sigs.astype('i2'))

mp.tight_layout()
mp.show()

　　

總結

以上是生活随笔為你收集整理的numpy之傅里叶定理的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。