一、基於傅里葉定理，用一組正弦函數合成方波
''' 三角函數通用函數 傅里葉定理：任何一個周期性曲線，無論多么跳躍或者不規則，都可以被解析成一組光滑的正弦函數的疊加 ---應用：合成方波(即不規則的方波由一組光滑的正弦函數疊加合成的) 如：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()
