2021/03/29

DSP 應用

數位音樂合成 digital music synthesis、數位語音合成 digital speech synthesis、數位語音辨識 digital speech recognition

數位音樂合成 digital music synthesis

使用訊號產生技術,包含週期性訊號及非週期性訊號,並透過數位訊號的排列組合,產生數位音樂

音樂的基本概念

音樂的基本構成要素:音高 pitch、節拍 beats、節奏 tempo

聲音的高低稱為音高 pitch,唱名:Do, Re, Mi, Fa, So, La, Si,對應的英語音名:C, D, E, F, G, A, B

鋼琴鍵盤的排列方式,是依照音高的順序,以中央 C 為基準向左右延伸,兩個相同音名鍵盤之間,有 8 個鍵盤,因此稱為八度音 octave。相鄰白鍵是相差一個全音,相鄰的白鍵與黑鍵,相差一個半音。

音高、音頻對照表:頻率,單位為赫茲。括號內為距離中央C(261.63赫茲)的半音距離。

0 1 2 3 4 5 6 7 8 9
C 16.352 (−48) 32.703 (−36) 65.406 (−24) 130.81 (−12) 261.63 (0) 523.25 (+12) 1046.5 (+24) 2093.0 (+36) 4186.0 (+48) 8372.0 (+60)
C♯/D♭ 17.324 (−47) 34.648 (−35) 69.296 (−23) 138.59 (−11) 277.18 (+1) 554.37 (+13) 1108.7 (+25) 2217.5 (+37) 4434.9 (+49) 8869.8 (+61)
D 18.354 (−46) 36.708 (−34) 73.416 (−22) 146.83 (−10) 293.66 (+2) 587.33 (+14) 1174.7 (+26) 2349.3 (+38) 4698.6 (+50) 9397.3 (+62)
D♯/E♭ 19.445 (−45) 38.891 (−33) 77.782 (−21) 155.56 (−9) 311.13 (+3) 622.25 (+15) 1244.5 (+27) 2489.0 (+39) 4978.0 (+51) 9956.1 (+63)
E 20.602 (−44) 41.203 (−32) 82.407 (−20) 164.81 (−8) 329.63 (+4) 659.26 (+16) 1318.5 (+28) 2637.0 (+40) 5274.0 (+52) 10548 (+64)
F 21.827 (−43) 43.654 (−31) 87.307 (−19) 174.61 (−7) 349.23 (+5) 698.46 (+17) 1396.9 (+29) 2793.8 (+41) 5587.7 (+53) 11175 (+65)
F♯/G♭ 23.125 (−42) 46.249 (−30) 92.499 (−18) 185.00 (−6) 369.99 (+6) 739.99 (+18) 1480.0 (+30) 2960.0 (+42) 5919.9 (+54) 11840 (+66)
G 24.500 (−41) 48.999 (−29) 97.999 (−17) 196.00 (−5) 392.00 (+7) 783.99 (+19) 1568.0 (+31) 3136.0 (+43) 6271.9 (+55) 12544 (+67)
G♯/A♭ 25.957 (−40) 51.913 (−28) 103.83 (−16) 207.65 (−4) 415.30 (+8) 830.61 (+20) 1661.2 (+32) 3322.4 (+44) 6644.9 (+56) 13290 (+68)
A 27.500 (−39) 55.000 (−27) 110.00 (−15) 220.00 (−3) 440.00 (+9) 880.00 (+21) 1760.0 (+33) 3520.0 (+45) 7040.0 (+57) 14080 (+69)
A♯/B♭ 29.135 (−38) 58.270 (−26) 116.54 (−14) 233.08 (−2) 466.16 (+10) 932.33 (+22) 1864.7 (+34) 3729.3 (+46) 7458.6 (+58) 14917 (+70)
B 30.868 (−37) 61.735 (−25) 123.47 (−13) 246.94 (−1) 493.88 (+11) 987.77 (+23) 1975.5 (+35) 3951.1 (+47) 7902.1 (+59) 15804 (+71)

樂曲中,每一個音都有自己的節拍 beats,代表這個音的時間長短,在五線譜中,節拍是用 音符 notes 表示,包含:全音符、二分音符、四分音符等等。ex: 五線譜中的拍號為 C 或 4/4,代表每一個小節有 4 拍,全音符代表這個音佔滿整個小節,因此為 4 拍,二分音符是全音符的一半,是 2 拍

另一個元素是節奏 tempo,就是音樂的快慢或速度。目前節奏通常是以每分鐘的節拍數 beats per minutes 決定。音樂的節奏包含:慢板、行板、中板、快板,與音樂要表達的情感有關。

小蜜蜂

import numpy as np
import wave
import struct

# 音符:音高 pitch + 節拍 beat
def note( pitch, beat ):
    fs = 44000
    amplitude = 30000
    # C, D, E, F, G, A, B 的頻率
    frequency = np.array( [ 261.6, 293.7, 329.6, 349.2, 392.0, 440.0, 493.9 ] )
    num_samples = beat * fs

    t = np.linspace( 0, beat, num_samples, endpoint = False )
    # 淡出效果
    a = np.linspace( 0, 1, num_samples, endpoint = False )

    # 弦波
    x = amplitude * a * np.cos( 2 * np.pi * frequency[ pitch - 1 ] * t )
    return x

def main():
    file = "little_bee.wav" # 檔案名稱

    # 音高 pitch
    pitches = np.array( [ 5, 3, 3, 4, 2, 2, 1, 2, 3, 4, 5, 5, 5,    \
                          5, 3, 3, 4, 2, 2, 1, 3, 5, 5, 3,          \
                          2, 2, 2, 2, 2, 3, 4, 3, 3, 3, 3, 3, 4, 5, \
                          5, 3, 3, 4, 2, 2, 1, 3, 5, 5, 1 ] )

    # 節拍 beat
    beats = np.array( [ 1, 1, 2, 1, 1, 2, 1, 1, 1, 1, 1, 1, 2,    \
                        1, 1, 2, 1, 1, 2, 1, 1, 1, 1, 4,          \
                        1, 1, 1, 1, 1, 1, 2, 1, 1, 1, 1, 1, 1, 2, \
                        1, 1, 2, 1, 1, 2, 1, 1, 1, 1, 4 ] )

    tempo = 0.5                 # 節奏(每拍0.5秒)
    fs = 44000

    # 時間總長度 = 節拍總和 * 節奏
    duration = sum( beats ) * tempo
    # 總樣本數 = 時間總長度 * 頻率
    num_samples = int( duration * fs )

    num_channels = 1            # 通道數
    samwidth = 2                # 樣本寬度
    num_frames = num_samples    # 音框數 = 樣本數
    comptype = "NONE"           # 壓縮型態
    compname = "not compressed" # 無壓縮

    num_notes = np.size( pitches )

    y = np.array( [ ] )
    for i in range( num_notes ):
        x = note( pitches[i], beats[i] * tempo )
        y = np.append( y, x )

    wav_file = wave.open( file, 'w' )
    wav_file.setparams(( num_channels, samwidth, fs, num_frames, comptype, compname ))

    for s in y:
        wav_file.writeframes( struct.pack( 'h', int( s ) ) )

    wav_file.close( )

main()

數位語音合成 digital speech synthesis

也就是 TTS, Text to Speech 的技術。

python 套件

  • Pyttsx Text to Speech
  • gTTS Text to Speech
  • eSpeak

數位語音辨識 digital speech recognition

Speech To Text

發展的技術:隱藏式馬可夫模型(Hidden Markov Models)、動態時間扭曲(Dynamic Time Warping, DTW)、人工神經網路(Artificial Neural Networks)、深度學習 (Deep Learning)、點對點自動語音辨識 (End-to-End Automatic Speech Recognition)。語音辨識率的準確度受到許多因素影響:雜訊、男生/女生、成人/兒童、口音、語意

新的人工智慧技術:遞迴神經網路 (Recurrent Neural Network, RNN),同時結合長短期記憶 (Long-Short Term Memory, LSTM)的技術最具代表性。AI 技術將成為語音辨識的主流

python 語音辨識 library: SpeechRecognition,支援許多 engines/apis

  • CMU Sphinx
  • Google Speech Recognition
  • Google Cloud Speech API
  • Wit.ai
  • Microsoft Bing Voice Recognition
  • Houndify API
  • IBM Speech to Txt
  • Snowboy Hotword Detection

References

數位訊號處理:Python程式實作(附範例光碟)(第二版)

2021/03/22

時頻分析

時間-頻率分析 Time-Frequency Analysis,首先介紹時頻分析的概念及數學工具 短時間傅立葉轉換 (Short-Time Fourier Transform, STFT),最後介紹範例與結果,用時頻圖 (spectrogram) 表示。

基本概念

通常輸入訊號是 non-stationary 訊號,會跟著時間改變。也就是頻率分量通常會隨著時間改變,因此無法在單一的頻譜分析非靜態(non-stationary)的訊號

例如一首歌曲,在不同時間點的頻率,會隨著時間改變,進而組成一首歌曲。人類發出的語音訊號,說話時每個字的音頻不同,進而組成一個句子。

為了分析或特性化 non-stationary 訊號,需要使用時間頻率分析技術,簡稱時頻分析技術

時頻分析 Time-Frequency Analysis:是訊號分析技術,同時包含時間域與頻率域,以時頻表示法Time-Frequency Representations 呈現。

因為時頻分析包含時間域與頻率域,可表示成二維訊號,稱為時頻表示法 Time-Frequency Representations ,通常是以二維圖形呈現

短時間傅立葉轉換

Short-Time Fourier Transform, STFT

STFT 定義:

\(STFT\{x(t)\}(τ, ω) = X(τ, ω) = \int_{-∞}^{∞}x(t)w(t-τ)e^{-jωt} {\rm d}t\)

其中 \(w(t)\) 為窗函數,通常是以原點為中心,包含 rectangular, hanning, gaussian window function 等等

注意 \(w, ω\) 代表不同的意義


離散短時間傅立葉轉換 Discrete STFT

定義:

\(STFT\{x[n]\}(m, ω) = X(m, ω) = \sum_{-∞}^{∞}x[n]w(n-m)e^{-jωn}\)

其中 \(w[n]\) 為窗函數

雖然 \(ω\) 是連續的,但因為 DSP 是用 快速傅立葉轉換 FFT 進行運算,因此頻率 \(ω\) 也會是離散的資料


時頻分析示意圖:

先將輸入訊號分成不同時間的區塊,通常是短暫且固定的時間,透過窗函數的運算,取得訊號區塊,每個訊號區塊,在經過傅立葉轉換,運算結果為複數,構成二維矩陣。最後,組合不同時間點訊號區塊的頻譜分析結果。

時頻圖 Spectrogram

定義:\(|X(m, ω)|^2\),其中 \(X(m, ω)\) 是 STFT 的結果

將 \(X(m, ω)\) 的二維複數矩陣結果,取強度 (magnitude) 平方,則形成二維的時數矩陣,就可以用二維的圖形呈現,稱為 spectrogram

STFT 主要缺點是解析度問題,由於窗函數的寬度固定,取得訊號區塊大小不同,使得時間域與頻率域的解析度也有所不同。

如訊號區塊的寬度較窄,則時間域解析度較佳,但在頻率域的解析度較差 (左圖)。

ex: 若數位訊號由弦波組成,時間長度為 10s,第一秒頻率 20Hz,第二秒頻率 40Hz,第三秒頻率 60Hz,遞增至 200Hz,取樣頻率為 1000Hz。假設窗函數為 rectangular window function,區塊長度為 100, 200, 500, 1000,求訊號的 STFT,以 spectrogram 表示。

import numpy as np
import scipy.signal as signal
import matplotlib.pyplot as plt

print( "Short-Time Fourier Transform" )
# n = eval( input( "Enter the length of segment: " ) )

def plot(x, fs, n, plotpos):
    f, t, Zxx = signal.stft( x, fs, window = 'boxcar', nperseg = n )

    plt.subplot(plotpos)
    plt.pcolormesh( t, f, abs(Zxx) )
    plt.xlabel( 'Time(Second)'+' (n='+str(n)+')' )
    plt.ylabel( 'Frequency(Hz)' )

fs = 1000
t = np.linspace( 0, 1, fs )

x = np.array( [ ] )
for i in range( 10 ):
    segment = np.cos( 2 * np.pi * ( ( i + 1 ) * 20 ) * t )
    x = np.append( x, segment )

# n: the length of segment
n = 100
plot(x, fs, n, '221')

n = 200
plot(x, fs, n, '222')

n = 500
plot(x, fs, n, '223')

n = 1000
plot(x, fs, n, '224')

plt.show( )

當弦波頻率增加,由於區塊長度不同,讓時間域與頻率域的解析度結果不同。當區塊長度增加,時間域的解析度越來越差,但頻率域的解析度越來越好。

時頻圖通常用 colormap 色彩圖呈現,亮度較高代表強度越大。

目前有一種 wavelet transform,可克服 STFT 的解析度問題。在高頻範圍,時間解析度較好,在低頻範圍,頻率解析度較好。


f, t, Zxx = signal.stft( x, fs, window = 'boxcar', nperseg = n ) 的部分可改用 SciPy 的 spectrogram 函式

f, t, Zxx = signal.spectrogram( x, fs )

import numpy as np
import scipy.signal as signal
import matplotlib.pyplot as plt

print( "Short-Time Fourier Transform" )
# n = eval( input( "Enter the length of segment: " ) )

def plot(x, fs):
    f, t, Zxx = signal.spectrogram( x, fs )

    plt.pcolormesh( t, f, abs(Zxx) )
    plt.xlabel( 'Time(Second)' )
    plt.ylabel( 'Frequency(Hz)' )

fs = 1000
t = np.linspace( 0, 1, fs )

x = np.array( [ ] )
for i in range( 10 ):
    segment = np.cos( 2 * np.pi * ( ( i + 1 ) * 20 ) * t )
    x = np.append( x, segment )

plot(x, fs)

plt.show( )

結果為

wav 的時頻分析

STFT 時頻圖

wav file 的 STFT 時頻圖

import numpy as np
import wave
from scipy.io import wavfile
import struct
import scipy.signal as signal
import matplotlib.pyplot as plt

infile  = input( "Input File: " )
fs, x = wavfile.read( infile )

# nperseg: 區塊長度 1000
f, t, Zxx = signal.stft( x, fs, nperseg = 1000 )

plt.pcolormesh( t, f, abs( Zxx ) )
plt.xlabel( 'Time(Second)' )
plt.ylabel( 'Frequency(Hz)' )

plt.show( )

sample rate: 11025Hz 的 wav file

SciPy 時頻圖(spectrogram)

import numpy as np
import wave
from scipy.io import wavfile
import struct
import scipy.signal as signal
import matplotlib.pyplot as plt

infile  = input( "Input File: " )   
fs, x = wavfile.read( infile )  
f, t, Zxx = signal.spectrogram( x, fs )

plt.pcolormesh( t, f, abs( Zxx ) )
plt.xlabel( 'Time(Second)' )
plt.ylabel( 'Frequency(Hz)' )

plt.show( )

跟上面的例子一樣的語音檔

References

數位訊號處理:Python程式實作(附範例光碟)(第二版)

2021/03/15

濾波器設計

首先介紹濾波器的基本概念:窗函數。濾波器的設計分成 FIR, IIR 兩種,其中 FIR Filter 使用 Window Method,IIR Filter 則包含 Butterworth Filter, Chebyshev Type I Filter, Chebyshev Type II Filter, 及 Elliptic Filter

濾波器在 DSP 系統中,可用來選取或抑制某些頻率範圍的輸入訊號,藉以產生符合需求的輸出訊號。

由於有限長度的脈衝響應無法實現 ideal filter,在設計實際的濾波器時,通常是考慮濾波器的頻率響應是否符合事先定義的規格,同時容許頻率響應有些微誤差。

實際的 lowpass filter 其頻率響應應如下圖,通過頻帶 passband 的範圍介於 \(0 ~ ω_p\) 之間,理想的頻率響應為 1;抑制頻帶 stopband 的範圍介於 \(ω_s ~ \pi\) 之間,理想的頻率響應為 0。其中 \(ω_p\) 稱為通過頻帶邊緣頻率 Passband Edge Frequency,\(ω_s\) 稱為抑制頻帶邊緣頻率 Stopband Edge Frequency。

實際的 lowpass filter 的設計,通常在通過頻帶或抑制頻帶內,容許頻率響應有些許波動(誤差),其中\(δ_p\) 稱為通過頻帶波紋 (Passband Ripple),\(δ_s\) 稱為抑制頻帶波紋 (Stopband Ripple)。在 \(ω_p ~ ω_s\) 之間,則是容許頻率響應從 1 逐漸降為 0,稱為過渡頻帶 (Transition Band)

filter 的設計,主要包含以下步驟:

  1. 定義濾波器規格:根據 DSP 應用的需求,選取濾波器的種類,例如:lowpass, highpass, bandpass, bandstop 等等。然後定義濾波器的規格,包含:passbadn 與 stopband 的截止頻率、取樣頻率等等參數

  2. 選取濾波器的種類:選取 FIR 或 IIR 濾波器作為設計目標。通常 FIR Filter 為穩定系統,若選取 IIR Filter ,則需要進一步檢驗其穩定性。

  3. 求轉換函式:根據濾波器規格,求濾波器的轉換函式 H(z),例如 FIR Filter 的轉換函式為

    \(H(z) = b_0+b_1z^{-1}+...+b_Mz^{-M}\)

    IIR 濾波器的轉換函式為:

    \(H(z) = \frac{b_0+b_1z^{-1}+...+b_Mz^{-M}}{a_0+a_1z^{-1}+...+a_Nz^{-N}}\)

    這個步驟的目的是求濾波器的係數:\(\{a_k\} k=0,1,...,N\) 與 \(\{b_k\} k=0,1,...,M\)

  4. 實現濾波器:在完成設計後,可採用硬體或軟體方式實現。

Window Function

idea lowpass filter 在時間域的脈衝響應為無限序列,因此,為了實現有限的脈衝響應,同時近似理想的 lowpass filter,最直接的方法是 Window Method。也就是定義 Window Function,用來擷取有限的脈衝響應。

Window Function 定義為:介於某特定區間的數學函示,區間外的數值為 0

任意函數與 window function 相乘後,結果就像是透過窗戶觀察該函數一樣。在 DSP,window function 經常用來進行頻譜分析、濾波器設計、音頻壓縮等等。

以下是幾個典型的 window function,這幾種 window function 都具有對稱性,SciPy 的Signal 有支援這些 functions。在濾波器設計過程中,窗函數通常是先經過平移,平移後以原點為中心,且對原點對稱,再與無限脈衝響應進行運算,藉以擷取有限的離散序列。

  • 矩形窗 Rectangular Window

    \( W(n) = \left\{ \begin{array}{ll} 1 & \mbox{if 0 ≤ n ≤ M-1} \\ 0 & \mbox{otherwise} \end{array} \right.\)

    其中 window 的大小為 M。由於矩形窗的形狀如同箱型車 boxcar,因此也被稱為 boxcar function

  • Hamming 窗

    \( W(n) = \left\{ \begin{array}{ll} 0.54-0.46cos(\frac{2 \pi n}{M-1}) & \mbox{if 0 ≤ n ≤ M-1} \\ 0 & \mbox{otherwise} \end{array} \right.\)

    這是根據 Richard W. Hamming 而命名的

  • Hanning 窗

    \( W(n) = \left\{ \begin{array}{ll} 0.5-0.5cos(\frac{2 \pi n}{M-1}) & \mbox{if 0 ≤ n ≤ M-1} \\ 0 & \mbox{otherwise} \end{array} \right.\)

    這是根據 Julus von Hann 而命名的,形狀跟 Hamming 窗相似

  • Bartlet 窗

    \( W(n) = \left\{ \begin{array}{ll} \frac{2}{M-1}(\frac{M-1}{2} - |n-\frac{M-1}{2}|) & \mbox{if 0 ≤ n ≤ M-1} \\ 0 & \mbox{otherwise} \end{array} \right.\)

    Bartlet window 的形狀是三角形

  • Blackman 窗

    \( W(n) = \left\{ \begin{array}{ll} 0.42-0.5cos(\frac{2 \pi n}{M-1})+0.08cos(\frac{4 \pi n}{M-1}) & \mbox{if 0 ≤ n ≤ M-1} \\ 0 & \mbox{otherwise} \end{array} \right.\)

    這是根據 Ralph B. Blackmand 而命名的

  • Kaiser 窗

    \(W(n) = I_0(β\sqrt{1-\frac{4n^2}{(M-1)^2}})/I_0(β), -\frac{M-1}{2} ≤ n ≤ \frac{M-1}{2}\)

    這是根據 Jim Kaiser 命名的。

    \(I_0\) 是修正的零階 Bessel 函式,參數 β 可用來近似其他窗函式。例如:β=0 為矩形窗、β=5 近似 Hamming 窗、β=6 近似 Hanning 窗、β=8.6 近似 Blackman 窗、β=14 則是 Kaiser 窗的建議初始值

統一選取 M=2*32+1 = 65,為有限脈衝響應的長度。window function 為對稱圖形,可平移到以原點為中心。

import numpy as np
import scipy.signal as signal
import matplotlib.pyplot as plt

M = 65
w1 = signal.boxcar( M )
w2 = signal.hamming( M )
w3 = signal.hann( M )
w4 = signal.bartlett( M )
w5 = signal.barthann( M )
w6 = signal.kaiser( M, 14 )

plt.subplot(231)
plt.plot( w1 )
plt.xlabel( 'n (rectangular, boxcar)' )
plt.ylabel( 'Amplitude' )

plt.subplot(232)
plt.plot( w2 )
plt.xlabel( 'n (hamming)' )
plt.ylabel( 'Amplitude' )

plt.subplot(233)
plt.plot( w3 )
plt.xlabel( 'n (hanning)' )
plt.ylabel( 'Amplitude' )

plt.subplot(234)
plt.plot( w4 )
plt.xlabel( 'n (bartlett)' )
plt.ylabel( 'Amplitude' )

plt.subplot(235)
plt.plot( w5 )
plt.xlabel( 'n (blackman)' )
plt.ylabel( 'Amplitude' )

plt.subplot(236)
plt.plot( w6 )
plt.xlabel( 'n (kaiser)' )
plt.ylabel( 'Amplitude' )

plt.show( )

FIR 濾波器設計

最簡單的方法是 window method,步驟如下

  1. 選取濾波器的種類:包含 lowpass, highpass, bandpass, bandstop 濾波器
  2. 定義濾波器的規格:如果是 lowpass, highpass,需定義截止頻率 (cutoff frequency),以 Hz 為單位,如果是 bandpass, bandstop,需定義兩個截止頻率
  3. 選取取樣頻率:sampling frequency,以 Hz 為單位。選取原則以輸入數位訊號為主,且對應的 Nyquist 頻率需大於上述的截止頻率
  4. 套用窗法:選取 window function。在此需選取濾波器的長度,也就是離散時間域的脈衝響應樣本數 (number of taps)
  5. FIR 濾波器係數:根據濾波器的種類,計算 FIR filter 的係數,並套用 window function,藉以擷取有限長度的脈衝響應。理論上,樣本數越多,越接近 ideal filter
  6. 頻率響應:根據脈衝響應,求其在頻率域的頻率響應,包含強度頻譜與相位頻譜,用來觀察濾波器的設計,是否符合原先地定義的規格需求

ex: 試用 window method 設計濾波器規格為

(1) lowpass filter:截止頻率 \(f_c=250Hz\)

(2) highpass filter:截止頻率 \(f_c=250Hz\)

(3) bandpass filter:截止頻率 \(f_1=200Hz, f_2=300Hz\)

(4) bandstop filter:截止頻率 \(f_1=200Hz, f_2=300Hz\)

並顯示強度頻率與相位頻譜


以下是用不同 window method 實作 lowpass filter 的結果

import numpy as np
import scipy.signal as signal
import matplotlib.pyplot as plt

def fir_filter(n, cutoff, win, pass_zero, freq, plotpos1, plotpos2):
    h = signal.firwin( n, cutoff, window = win, pass_zero = pass_zero, fs = freq )

    w, H = signal.freqz( h )
    magnitude = abs( H )
    phase = np.angle( H )

    plt.subplot(plotpos1)
    plt.plot( w, magnitude )
    plt.xlabel( r'$\omega$ '+'強度頻譜 ('+str(win)+')' )
    plt.ylabel( 'Magnitude' )

    plt.subplot(plotpos2)
    plt.plot( w, phase )
    plt.xlabel( r'$\omega$ '+'相位頻譜 ('+str(win)+')' )
    plt.ylabel( 'Phase' )

# cutoff frequency(Hz)
cutoff = 250
# numeber of taps
n = 31
# pass_zero = {True, False, 'bandpass', 'lowpass', 'highpass', 'bandstop'}
pass_zero = 'lowpass'
# sampling frequency (Hz)
freq = 1000


# 替換中文字型,處理 matplotlib 中文 label 問題
# https://blog.csdn.net/gmr2453929471/article/details/78655834
plt.rcParams['font.sans-serif'] = ['Heiti TC'] # 步驟一(替換sans-serif字型)
plt.rcParams['axes.unicode_minus'] = False  # 步驟二(解決座標軸負數的負號顯示問題)

plt.figure( 1 )

# rectangular (boxcar) window method
win='boxcar'
fir_filter(n, cutoff, win, pass_zero, freq, 231, 234)

# hamming window method
win='hamming'
fir_filter(n, cutoff, win, pass_zero, freq, 232, 235)

# hanning window method
win='hanning'
fir_filter(n, cutoff, win, pass_zero, freq, 233, 236)

plt.figure( 2 )

# bartlett window method
win='bartlett'
fir_filter(n, cutoff, win, pass_zero, freq, 231, 234)

# blackman window method
win='blackman'
fir_filter(n, cutoff, win, pass_zero, freq, 232, 235)

# kaiser window method
win=('kaiser',14)
fir_filter(n, cutoff, win, pass_zero, freq, 233, 236)

plt.show( )

使用 rectangular window method 設計的 lowpass filter,其頻率響應有 Gibbs 現象,結果較不理想。使用 blackman 或kaiser 設計的 lowpass filter,結果較理想。

截止角頻率經過正規化為 \(2 \pi f_c/f_s = 2 \pi (250) / 1000 = \pi/2\),角頻率介於 \(0~\pi\) 之間


以下是用 hamming window method 實作四種 filter 的結果

import numpy as np
import scipy.signal as signal
import matplotlib.pyplot as plt

def fir_filter(n, cutoff, win, pass_zero, freq, plotpos1):
    h = signal.firwin( n, cutoff, window = win, pass_zero = pass_zero, fs = freq )

    w, H = signal.freqz( h )
    magnitude = abs( H )
    phase = np.angle( H )

    plt.subplot(plotpos1)
    plt.plot( w, magnitude )
    plt.xlabel( r'$\omega$ '+'('+str(pass_zero)+')' )
    plt.ylabel( 'Magnitude' )

    # plt.subplot(plotpos2)
    # plt.plot( w, phase )
    # plt.xlabel( r'$\omega$ '+'相位頻譜 ('+str(win)+')' )
    # plt.ylabel( 'Phase' )

# cutoff frequency(Hz)
cutoff = 250
f1 = 200
f2 = 300

# numeber of taps
n = 31
# sampling frequency (Hz)
freq = 1000

# hamming window method
win='hamming'

plt.figure( 1 )

# lowpass
# pass_zero = {True, False, 'bandpass', 'lowpass', 'highpass', 'bandstop'}
pass_zero='lowpass'
fir_filter(n, cutoff, win, pass_zero, freq, 221)

pass_zero='highpass'
fir_filter(n, cutoff, win, pass_zero, freq, 222)

pass_zero='bandpass'
fir_filter(n, [f1, f2], win, pass_zero, freq, 223)

pass_zero='bandstop'
fir_filter(n, [f1, f2], win, pass_zero, freq, 224)

plt.show( )

IIR 濾波器設計

IIR filters 分為以下幾種

  1. Butterworth Filter
  2. Chebyshev Type I Filter
  3. Chebyshev Type II Filter
  4. Elliptic Filter 橢圓濾波器

IIR filter 的設計步驟,跟 FIR filter 的步驟類似

  1. 選取濾波器的種類:包含 lowpass, highpass, bandpass, bandstop 濾波器

  2. 定義濾波器的規格:如果是 lowpass, highpass,需定義截止頻率 (cutoff frequency),以 Hz 為單位,如果是 bandpass, bandstop,需定義兩個截止頻率。另外需要定義通過頻帶波紋與抑制頻帶波紋,單位為 dB

  3. 選取取樣頻率:sampling frequency,以 Hz 為單位。選取原則以輸入數位訊號為主,且對應的 Nyquist 頻率需大於上述的截止頻率

  4. 濾波器參數:根據濾波器的規格,計算濾波器的階數 (Order) 與 -3dB 截止頻率

  5. IIR 濾波器係數:根據濾波器的種類及相關參數,計算轉換函式

    \(H(z)=\frac{Y(z)}{X(z)} = \frac{b_0+b_1z^{-1}+...+b_Mz^{-M}}{a_0+z_1z^{-1}+...+a_Nz^{-N}}\)

    也就是計算濾波器的係數 \(\{a_k\}, k=1,2,...,N\) \(\{b_k\}, k=0,1,...M\)

  6. 頻率響應:根據脈衝響應,求其在頻率域的頻率響應,包含強度頻譜與相位頻譜,用來觀察濾波器的設計,是否符合原先地定義的規格需求

Butterworth 濾波器

是一種通過頻帶的頻率響應非常平坦的 DSP,由英國工程師 Stephen Butterworth 提出

根據下列規格,設計 Butterworth filters

(1) lowpass filter & highpass filter

通過頻帶邊緣頻率 200 Hz

抑制頻帶邊緣頻率 300 Hz

通過頻帶波紋 0.5 dB

抑制頻帶波紋 50 dB

取樣頻率 1000Hz

(2) bandpass & bandstop filter

通過頻帶邊緣頻率 200, 300 Hz

抑制頻帶邊緣頻率 100, 400 Hz

通過頻帶波紋 0.5 dB

抑制頻帶波紋 50 dB

取樣頻率 1000Hz

import numpy as np
import scipy.signal as signal
import matplotlib.pyplot as plt

def plot(a, b, plotname, plotpos1, plotpos2):
    w, H = signal.freqz( b, a )
    magnitude = abs( H )
    phase = np.angle( H )

    plt.subplot(plotpos1)
    plt.plot( w, magnitude )
    plt.xlabel( r'$\omega$ 強度頻譜('+plotname+')' )
    plt.ylabel( 'Magnitude' )

    plt.subplot(plotpos2)
    plt.plot( w, phase )
    plt.xlabel( r'$\omega$ 相位頻譜('+plotname+')' )
    plt.ylabel( 'Phase' )

# 替換中文字型,處理 matplotlib 中文 label 問題
# https://blog.csdn.net/gmr2453929471/article/details/78655834
plt.rcParams['font.sans-serif'] = ['Heiti TC'] # 步驟一(替換sans-serif字型)
plt.rcParams['axes.unicode_minus'] = False  # 步驟二(解決座標軸負數的負號顯示問題)


# lowpass filter
# passband edge frequency(Hz)
fp = 200
# stopband edge frequency(Hz)
fs = 300
# passband ripple(dB)
rp = 0.5
# stopband ripple(dB)
rs = 50
# sampling frequency
Fs = 1000

wp = 2 * fp / Fs
ws = 2 * fs / Fs

n, wn = signal.buttord( wp, ws, rp, rs )
b, a = signal.butter( n, wn, 'lowpass' )

plt.figure( 1 )
plot(a, b, 'lowpass', '121', '122')


# highpass filter
# passband edge frequency(Hz)
fp = 200
# stopband edge frequency(Hz)
fs = 300
# passband ripple(dB)
rp = 0.5
# stopband ripple(dB)
rs = 50
# sampling frequency
Fs = 1000

wp = 2 * fp / Fs
ws = 2 * fs / Fs

n, wn = signal.buttord( wp, ws, rp, rs )
b, a = signal.butter( n, wn, 'highpass' )

plt.figure( 2 )
plot(a, b, 'highpass', '121', '122')


# bandpass filter
# 1st passband edge frequency(Hz)
fp1 = 200
# 2nd passband edge frequency(Hz)
fp2 = 300
# 1st stopband edge frequency(Hz)
fs1 = 100
# 2nd stopband edge frequency(Hz)
fs2 = 400
# passband ripple(dB)
rp = 0.5
# stopband ripple(dB)
rs = 50
# sampling frequency
Fs = 1000

wp1 = 2 * fp1 / Fs
wp2 = 2 * fp2 / Fs
ws1 = 2 * fs1 / Fs
ws2 = 2 * fs2 / Fs

n, wn = signal.buttord( [wp1, wp2], [ws1, ws2], rp, rs )
b, a = signal.butter( n, wn, 'bandpass' )

plt.figure( 3 )
plot(a, b, 'bandpass', '121', '122')


# bandstop filter
# 1st passband edge frequency(Hz)
fp1 = 200
# 2nd passband edge frequency(Hz)
fp2 = 300
# 1st stopband edge frequency(Hz)
fs1 = 100
# 2nd stopband edge frequency(Hz)
fs2 = 400
# passband ripple(dB)
rp = 0.5
# stopband ripple(dB)
rs = 50
# sampling frequency
Fs = 1000

wp1 = 2 * fp1 / Fs
wp2 = 2 * fp2 / Fs
ws1 = 2 * fs1 / Fs
ws2 = 2 * fs2 / Fs

n, wn = signal.buttord( [wp1, wp2], [ws1, ws2], rp, rs )
b, a = signal.butter( n, wn, 'bandstop' )

plt.figure( 4 )
plot(a, b, 'bandstop', '121', '122')

plt.show( )

Chebyshev Type I Filter

Chebyshev Filter 是一種在通過頻帶或抑制頻帶上頻率響應等波動的濾波器。在通過頻帶波動的濾波器為 Chebyshev Type I Filter,在抑制頻帶上波動的濾波器為 Chebyshev Type II Filter

根據下列規格,設計 Chebyshev Type I Filters

(1) lowpass filter & highpass filter

通過頻帶邊緣頻率 200 Hz

抑制頻帶邊緣頻率 300 Hz

通過頻帶波紋 0.5 dB

抑制頻帶波紋 50 dB

取樣頻率 1000Hz

(2) bandpass & bandstop filter

通過頻帶邊緣頻率 200, 300 Hz

抑制頻帶邊緣頻率 100, 400 Hz

通過頻帶波紋 0.5 dB

抑制頻帶波紋 50 dB

取樣頻率 1000Hz

import numpy as np
import scipy.signal as signal
import matplotlib.pyplot as plt

def plot(a, b, plotname, plotpos1, plotpos2):
    w, H = signal.freqz( b, a )
    magnitude = abs( H )
    phase = np.angle( H )

    plt.subplot(plotpos1)
    plt.plot( w, magnitude )
    plt.xlabel( r'$\omega$ 強度頻譜('+plotname+')' )
    plt.ylabel( 'Magnitude' )

    plt.subplot(plotpos2)
    plt.plot( w, phase )
    plt.xlabel( r'$\omega$ 相位頻譜('+plotname+')' )
    plt.ylabel( 'Phase' )

# 替換中文字型,處理 matplotlib 中文 label 問題
# https://blog.csdn.net/gmr2453929471/article/details/78655834
plt.rcParams['font.sans-serif'] = ['Heiti TC'] # 步驟一(替換sans-serif字型)
plt.rcParams['axes.unicode_minus'] = False  # 步驟二(解決座標軸負數的負號顯示問題)

# lowpass filter
# passband edge frequency(Hz)
fp = 200
# stopband edge frequency(Hz)
fs = 300
# passband ripple(dB)
rp = 0.5
# stopband ripple(dB)
rs = 50
# sampling frequency
Fs = 1000

wp = 2 * fp / Fs
ws = 2 * fs / Fs

n, wn = signal.cheb1ord( wp, ws, rp, rs )
b, a = signal.cheby1( n, rp, wn, 'lowpass' )

plt.figure( 1 )
plot(a, b, 'lowpass', '221', '222')


# highpass filter
# passband edge frequency(Hz)
fp = 200
# stopband edge frequency(Hz)
fs = 300
# passband ripple(dB)
rp = 0.5
# stopband ripple(dB)
rs = 50
# sampling frequency
Fs = 1000

wp = 2 * fp / Fs
ws = 2 * fs / Fs

n, wn = signal.cheb1ord( wp, ws, rp, rs )
b, a = signal.cheby1( n, rp, wn, 'highpass' )

plot(a, b, 'highpass', '223', '224')


# bandpass filter
# 1st passband edge frequency(Hz)
fp1 = 200
# 2nd passband edge frequency(Hz)
fp2 = 300
# 1st stopband edge frequency(Hz)
fs1 = 100
# 2nd stopband edge frequency(Hz)
fs2 = 400
# passband ripple(dB)
rp = 0.5
# stopband ripple(dB)
rs = 50
# sampling frequency
Fs = 1000

wp1 = 2 * fp1 / Fs
wp2 = 2 * fp2 / Fs
ws1 = 2 * fs1 / Fs
ws2 = 2 * fs2 / Fs

n, wn = signal.cheb1ord( [ wp1, wp2 ], [ ws1, ws2 ], rp, rs )
b, a = signal.cheby1( n, rp, wn, 'bandpass' )

plt.figure( 2 )
plot(a, b, 'bandpass', '221', '222')


# bandstop filter
# 1st passband edge frequency(Hz)
fp1 = 200
# 2nd passband edge frequency(Hz)
fp2 = 300
# 1st stopband edge frequency(Hz)
fs1 = 100
# 2nd stopband edge frequency(Hz)
fs2 = 400
# passband ripple(dB)
rp = 0.5
# stopband ripple(dB)
rs = 50
# sampling frequency
Fs = 1000

wp1 = 2 * fp1 / Fs
wp2 = 2 * fp2 / Fs
ws1 = 2 * fs1 / Fs
ws2 = 2 * fs2 / Fs

n, wn = signal.cheb1ord( [ wp1, wp2 ], [ ws1, ws2 ], rp, rs )
b, a = signal.cheby1( n, rp, wn, 'bandstop' )

plot(a, b, 'bandstop', '223', '224')

plt.show( )

Chebyshev Type II Filter

在抑制頻帶上波動的濾波器為 Chebyshev Type II Filter

根據下列規格,設計 Chebyshev Type II Filters

(1) lowpass filter & highpass filter

通過頻帶邊緣頻率 200 Hz

抑制頻帶邊緣頻率 300 Hz

通過頻帶波紋 0.5 dB

抑制頻帶波紋 50 dB

取樣頻率 1000Hz

(2) bandpass & bandstop filter

通過頻帶邊緣頻率 200, 300 Hz

抑制頻帶邊緣頻率 100, 400 Hz

通過頻帶波紋 0.5 dB

抑制頻帶波紋 50 dB

取樣頻率 1000Hz

import numpy as np
import scipy.signal as signal
import matplotlib.pyplot as plt

def plot(a, b, plotname, plotpos1, plotpos2):
    w, H = signal.freqz( b, a )
    magnitude = abs( H )
    phase = np.angle( H )

    plt.subplot(plotpos1)
    plt.plot( w, magnitude )
    plt.xlabel( r'$\omega$ 強度頻譜('+plotname+')' )
    plt.ylabel( 'Magnitude' )

    plt.subplot(plotpos2)
    plt.plot( w, phase )
    plt.xlabel( r'$\omega$ 相位頻譜('+plotname+')' )
    plt.ylabel( 'Phase' )

# 替換中文字型,處理 matplotlib 中文 label 問題
# https://blog.csdn.net/gmr2453929471/article/details/78655834
plt.rcParams['font.sans-serif'] = ['Heiti TC'] # 步驟一(替換sans-serif字型)
plt.rcParams['axes.unicode_minus'] = False  # 步驟二(解決座標軸負數的負號顯示問題)

# lowpass filter
# passband edge frequency(Hz)
fp = 200
# stopband edge frequency(Hz)
fs = 300
# passband ripple(dB)
rp = 0.5
# stopband ripple(dB)
rs = 50
# sampling frequency
Fs = 1000

wp = 2 * fp / Fs
ws = 2 * fs / Fs

n, wn = signal.cheb2ord( wp, ws, rp, rs )
b, a = signal.cheby2( n, rp, wn, 'lowpass' )

plt.figure( 1 )
plot(a, b, 'lowpass', '221', '222')


# highpass filter
# passband edge frequency(Hz)
fp = 200
# stopband edge frequency(Hz)
fs = 300
# passband ripple(dB)
rp = 0.5
# stopband ripple(dB)
rs = 50
# sampling frequency
Fs = 1000

wp = 2 * fp / Fs
ws = 2 * fs / Fs

n, wn = signal.cheb2ord( wp, ws, rp, rs )
b, a = signal.cheby2( n, rp, wn, 'highpass' )

plot(a, b, 'highpass', '223', '224')


# bandpass filter
# 1st passband edge frequency(Hz)
fp1 = 200
# 2nd passband edge frequency(Hz)
fp2 = 300
# 1st stopband edge frequency(Hz)
fs1 = 100
# 2nd stopband edge frequency(Hz)
fs2 = 400
# passband ripple(dB)
rp = 0.5
# stopband ripple(dB)
rs = 50
# sampling frequency
Fs = 1000

wp1 = 2 * fp1 / Fs
wp2 = 2 * fp2 / Fs
ws1 = 2 * fs1 / Fs
ws2 = 2 * fs2 / Fs

n, wn = signal.cheb2ord( [ wp1, wp2 ], [ ws1, ws2 ], rp, rs )
b, a = signal.cheby2( n, rp, wn, 'bandpass' )

plt.figure( 2 )
plot(a, b, 'bandpass', '221', '222')


# bandstop filter
# 1st passband edge frequency(Hz)
fp1 = 200
# 2nd passband edge frequency(Hz)
fp2 = 300
# 1st stopband edge frequency(Hz)
fs1 = 100
# 2nd stopband edge frequency(Hz)
fs2 = 400
# passband ripple(dB)
rp = 0.5
# stopband ripple(dB)
rs = 50
# sampling frequency
Fs = 1000

wp1 = 2 * fp1 / Fs
wp2 = 2 * fp2 / Fs
ws1 = 2 * fs1 / Fs
ws2 = 2 * fs2 / Fs

n, wn = signal.cheb2ord( [ wp1, wp2 ], [ ws1, ws2 ], rp, rs )
b, a = signal.cheby2( n, rp, wn, 'bandstop' )

plot(a, b, 'bandstop', '223', '224')

plt.show( )

Elliptic Filter

是一種在通過頻帶與抑制頻帶波紋的濾波器。橢圓濾波器比其他類型的濾波器,在階數相同的條件下,有比較小的通過頻帶與抑制頻帶波動。

根據下列規格,設計 Elliptic Filters

(1) lowpass filter & highpass filter

通過頻帶邊緣頻率 200 Hz

抑制頻帶邊緣頻率 300 Hz

通過頻帶波紋 0.5 dB

抑制頻帶波紋 50 dB

取樣頻率 1000Hz

(2) bandpass & bandstop filter

通過頻帶邊緣頻率 200, 300 Hz

抑制頻帶邊緣頻率 100, 400 Hz

通過頻帶波紋 0.5 dB

抑制頻帶波紋 50 dB

取樣頻率 1000Hz

import numpy as np
import scipy.signal as signal
import matplotlib.pyplot as plt

def plot(a, b, plotname, plotpos1, plotpos2):
    w, H = signal.freqz( b, a )
    magnitude = abs( H )
    phase = np.angle( H )

    plt.subplot(plotpos1)
    plt.plot( w, magnitude )
    plt.xlabel( r'$\omega$ 強度頻譜('+plotname+')' )
    plt.ylabel( 'Magnitude' )

    plt.subplot(plotpos2)
    plt.plot( w, phase )
    plt.xlabel( r'$\omega$ 相位頻譜('+plotname+')' )
    plt.ylabel( 'Phase' )

# 替換中文字型,處理 matplotlib 中文 label 問題
# https://blog.csdn.net/gmr2453929471/article/details/78655834
plt.rcParams['font.sans-serif'] = ['Heiti TC'] # 步驟一(替換sans-serif字型)
plt.rcParams['axes.unicode_minus'] = False  # 步驟二(解決座標軸負數的負號顯示問題)

# lowpass filter
# passband edge frequency(Hz)
fp = 200
# stopband edge frequency(Hz)
fs = 300
# passband ripple(dB)
rp = 0.5
# stopband ripple(dB)
rs = 50
# sampling frequency
Fs = 1000

wp = 2 * fp / Fs
ws = 2 * fs / Fs

n, wn = signal.ellipord( wp, ws, rp, rs )
b, a = signal.ellip( n, rp, rs, wn, 'lowpass' )

plt.figure( 1 )
plot(a, b, 'lowpass', '221', '222')


# highpass filter
# passband edge frequency(Hz)
fp = 200
# stopband edge frequency(Hz)
fs = 300
# passband ripple(dB)
rp = 0.5
# stopband ripple(dB)
rs = 50
# sampling frequency
Fs = 1000

wp = 2 * fp / Fs
ws = 2 * fs / Fs

n, wn = signal.ellipord( wp, ws, rp, rs )
b, a = signal.ellip( n, rp, rs, wn, 'highpass' )

plot(a, b, 'highpass', '223', '224')


# bandpass filter
# 1st passband edge frequency(Hz)
fp1 = 200
# 2nd passband edge frequency(Hz)
fp2 = 300
# 1st stopband edge frequency(Hz)
fs1 = 100
# 2nd stopband edge frequency(Hz)
fs2 = 400
# passband ripple(dB)
rp = 0.5
# stopband ripple(dB)
rs = 50
# sampling frequency
Fs = 1000

wp1 = 2 * fp1 / Fs
wp2 = 2 * fp2 / Fs
ws1 = 2 * fs1 / Fs
ws2 = 2 * fs2 / Fs

n, wn = signal.ellipord( [ wp1, wp2 ], [ ws1, ws2 ], rp, rs )
b, a = signal.ellip( n, rp, rs, wn, 'bandpass' )

plt.figure( 2 )
plot(a, b, 'bandpass', '221', '222')


# bandstop filter
# 1st passband edge frequency(Hz)
fp1 = 200
# 2nd passband edge frequency(Hz)
fp2 = 300
# 1st stopband edge frequency(Hz)
fs1 = 100
# 2nd stopband edge frequency(Hz)
fs2 = 400
# passband ripple(dB)
rp = 0.5
# stopband ripple(dB)
rs = 50
# sampling frequency
Fs = 1000

wp1 = 2 * fp1 / Fs
wp2 = 2 * fp2 / Fs
ws1 = 2 * fs1 / Fs
ws2 = 2 * fs2 / Fs

n, wn = signal.ellipord( [ wp1, wp2 ], [ ws1, ws2 ], rp, rs )
b, a = signal.ellip( n, rp, rs, wn, 'bandstop' )

plot(a, b, 'bandstop', '223', '224')

plt.show( )

References

數位訊號處理:Python程式實作(附範例光碟)(第二版)

2021/03/08

頻率域 DSP

頻率域的數位訊號處理 Digital Signal Processing in Frequency Domain

到目前為止,DSP 是以離散時間域的數學運算為主軸,例如:卷積運算、移動平均濾波器、FIR 濾波器、IIR 濾波器。本章是以頻率域的數位訊號處理為主軸

頻率域 DSP 的系統方塊圖

處理步驟如下:

  1. 輸入的數位訊號 x[n],取離散時間傅立葉轉換 DTFT,結果為

    \(X(e^{jω}) = F\{x[n]\}\)

    在DSP 實務中,則是採用離散傅立葉轉換 (Discrete Fourier Transform, DFT),即

    \(X[k] = \sum_{n=0}^{N-1}x[n]e^{-j2 \pi kn/N}, k=0,1,...,N-1\)

    為了縮短處理時間,通常採用快速傅立葉轉換 Fast Fourier Transform, FFT 的演算法進行運算。注意運算結果為複數陣列

  2. 根據系統輸入/輸出的目的或規格,進行濾波器的設計,藉以產生頻率響應,即

    \(H(e^{jω}) = F\{h[n]\}\)

    在 DSP 實務中,根據濾波器的種類,設計濾波器的頻率響應 Frequency Response,通常頻率響應在特定頻率範圍的響應值(或強度值),介於 0~1 之間

  3. 套用濾波器,通常使用卷積定理,藉以產生輸出訊號的頻率域表示法

    \(Y(e^{jω}) = H(e^{jω})X(e^{jω})\)

    在 DSP 實務中,運用陣列的點對點乘法 (Pointwise Multiplication) 運算,達到濾波的效果,即

    \(Y[k] = H[k] \cdot X[k]\)

    結果通常是複數陣列

  4. 最後,取反離散時間傅立葉轉換(Inverse DTFT),即可得到輸出訊號

    \(y[n] = F^{-1}\{Y(e^{jω})\}\)

    在 DSP 實務中,採用反離散傅立葉轉換 Inverse DFT,即

    \(y[n] = \frac{1}{N} \sum_{k=0}^{N-1} Y[k]e^{j2 \pi kn/N}, n=0,1,2,...,N-1\)

    在此,也是使用反快速傅立葉轉換 Inverse FFT,縮短處理時間。由於處理後的結果為複數,通常是取實數部分(忽略虛數),作為輸出訊號,必要時,進一步將數值作正規化處理,藉以控制輸出訊號的數值範圍

理想濾波器

理想濾波器在離散時間域的脈衝響應是有限數列,因此無法在離散時間域實現理想濾波器,然而透過頻率域 DSP,可實現理想濾波器。

ex: 輸入訊號諧波定義如下

\(x(t) = cos(2 \pi \cdot 10 \cdot t)+ cos(2 \pi \cdot 20 \cdot t) + cos(2 \pi \cdot 30 \cdot t)\)

取樣頻率 \(f_s=500Hz\),試套用以下理想濾波器

(1) 理想低通濾波器,截止頻率 \(f_c=15Hz\)

(2) 理想高通濾波器,截止頻率 \(f_c=15Hz\)

(3) 理想帶通濾波器,截止頻率 \(f_1=15Hz, f_2=25Hz\)

(4) 理想帶阻濾波器,截止頻率 \(f_1=15Hz, f_2=25Hz\)

import numpy as np
from numpy.fft import fft, fftshift, ifft, fftfreq
import matplotlib.pyplot as plt

def ideal_lowpass_filtering( x, cutoff, fs ):
    X = fft( x )
    H = np.zeros( fs )
    for i in range( -cutoff, cutoff + 1 ):
        H[i] = 1
    Y = H * X
    y = ifft( Y )
    y = y.real
    return y

def ideal_highpass_filtering( x, cutoff, fs ):
    X = fft( x )
    H = np.zeros( fs )
    for i in range( -cutoff, cutoff + 1 ):
        H[i] = 1
    H = 1 - H
    Y = H * X
    y = ifft( Y )
    y = y.real
    return y

def ideal_bandpass_filtering( x, f1, f2, fs ):
    X = fft( x )
    H = np.zeros( fs )
    for i in range( f1, f2 + 1 ):
        H[i] = 1
    for i in range( -f1, -f2 - 1, -1 ):
        H[i] = 1
    Y = H * X
    y = ifft( Y )
    y = y.real
    return y

def ideal_bandstop_filtering( x, f1, f2, fs ):
    X = fft( x )
    H = np.zeros( fs )
    for i in range( f1, f2 + 1 ):
        H[i] = 1
    for i in range( -f1, -f2 - 1, -1 ):
        H[i] = 1
    H = 1 - H
    Y = H * X
    y = ifft( Y )
    y = y.real
    return y

def ideal_allpass_filtering( x ):
    X = fft( x )
    Y = X
    y = ifft( Y )
    y = y.real
    return y

def main( ):
    # print( "DSP in Frequency Domain" )
    # print( "(1) Ideal Lowpass Filtering" )
    # print( "(2) Ideal Highpass Filtering" )
    # print( "(3) Ideal Bandpass Filtering" )
    # print( "(4) Ideal Bandstop Filtering" )
    # print( "(5) Ideal Allpass Filtering" )

    # choice = eval( input( "Please enter your choice: " ) )

    # if choice == 1 or choice == 2:
    #   fc = eval( input( "Please enter cutoff frequency(Hz): " ) )

    # if choice == 3 or choice == 4:
    #   f1 = eval( input( "Please enter frequency f1(Hz): " ) )
    #   f2 = eval( input( "Please enter frequency f2(Hz): " ) )

    fc = 15

    f1 = 15
    f2 = 25

    # 取樣頻率
    fs = 500
    t = np.linspace( 0, 1, fs, endpoint = False )
    # 原始訊號 (諧波)
    x = np.cos( 2 * np.pi * 10 * t ) + np.cos( 2 * np.pi * 20 * t ) + np.cos( 2 * np.pi * 30 * t )

    # fftfreq: 離散傅立葉轉換的採樣頻率
    # fftshift: 將 0 移到正中心
    f = fftshift( fftfreq( fs, 1 / fs ) )
    Xm = abs( fftshift( fft( x ) ) )

    # 替換中文字型,處理 matplotlib 中文 label 問題
    # https://blog.csdn.net/gmr2453929471/article/details/78655834
    plt.rcParams['font.sans-serif'] = ['Heiti TC'] # 步驟一(替換sans-serif字型)
    plt.rcParams['axes.unicode_minus'] = False  # 步驟二(解決座標軸負數的負號顯示問題)

    # 輸入訊號圖形
    plt.figure( 1 )
    plt.subplot(121)
    plt.plot( x )
    plt.xlabel( 't (second) (輸入訊號)' )
    plt.ylabel( 'Amplitude' )

    plt.subplot(122)
    plt.plot( f, Xm )
    plt.xlabel( 'f (輸入訊號頻譜)' )
    plt.ylabel( 'Magnitude' )

    # if choice == 1:
    #   y = ideal_lowpass_filtering( x, fc, fs )
    # elif choice == 2:
    #   y = ideal_highpass_filtering( x, fc, fs )
    # elif choice == 3:
    #   y = ideal_bandpass_filtering( x, f1, f2, fs )
    # elif choice == 4:
    #   y = ideal_bandstop_filtering( x, f1, f2, fs )
    # else:
    #   y = ideal_allpass_filtering( x )

    # ideal_lowpass_filtering
    y = ideal_lowpass_filtering( x, fc, fs )
    Ym = abs( fftshift( fft( y ) ) )

    plt.figure( 2 )
    plt.subplot(221)
    plt.plot( y )
    plt.xlabel( 't (second) (ideal_lowpass 輸出訊號)' )
    plt.ylabel( 'Amplitude' )

    plt.subplot(222)
    plt.plot( f, Ym )
    plt.xlabel( 'f  (ideal_lowpass 輸出訊號頻譜)' )
    plt.ylabel( 'Magnitude' )

    # ideal_highpass_filtering
    y = ideal_highpass_filtering( x, fc, fs )
    Ym = abs( fftshift( fft( y ) ) )

    plt.figure( 2 )
    plt.subplot(223)
    plt.plot( y )
    plt.xlabel( 't (second) (ideal_highpass 輸出訊號)' )
    plt.ylabel( 'Amplitude' )

    plt.subplot(224)
    plt.plot( f, Ym )
    plt.xlabel( 'f  (ideal_highpass 輸出訊號頻譜)' )
    plt.ylabel( 'Magnitude' )

    # ideal_bandpass_filtering
    y = ideal_bandpass_filtering( x, f1, f2, fs )
    Ym = abs( fftshift( fft( y ) ) )
    plt.figure( 3 )
    plt.subplot(221)
    plt.plot( y )
    plt.xlabel( 't (second) (ideal_bandpass 輸出訊號)' )
    plt.ylabel( 'Amplitude' )

    plt.subplot(222)
    plt.plot( f, Ym )
    plt.xlabel( 'f  (ideal_bandpass 輸出訊號頻譜)' )
    plt.ylabel( 'Magnitude' )

    # ideal_bandstop_filtering
    y = ideal_bandstop_filtering( x, f1, f2, fs )
    Ym = abs( fftshift( fft( y ) ) )
    plt.subplot(223)
    plt.plot( y )
    plt.xlabel( 't (second) (ideal_bandstop 輸出訊號)' )
    plt.ylabel( 'Amplitude' )

    plt.subplot(224)
    plt.plot( f, Ym )
    plt.xlabel( 'f  (ideal_bandstop 輸出訊號頻譜)' )
    plt.ylabel( 'Magnitude' )

    # ideal_allpass_filtering
    y = ideal_allpass_filtering( x )
    Ym = abs( fftshift( fft( y ) ) )
    plt.figure( 4 )
    plt.subplot(221)
    plt.plot( y )
    plt.xlabel( 't (second) (ideal_allpass 輸出訊號)' )
    plt.ylabel( 'Amplitude' )

    plt.subplot(222)
    plt.plot( f, Ym )
    plt.xlabel( 'f  (ideal_allpass 輸出訊號頻譜)' )
    plt.ylabel( 'Magnitude' )

    plt.show( )

main( )

輸入訊號有 10, 20, 30 Hz 三種頻率

ideal lowpass filter:抑制超過 15Hz 的頻率,輸出只剩下 10Hz 的頻率分量

ideal highpass filter:抑制 15Hz 以下的頻率,輸出 20, 30 Hz 的頻率分量

ideal bandpass filter:允許 15~25 Hz 的頻率分量,剩下 20 Hz

ideal bandstop filter:阻絕 15~25 Hz 的頻率分量,剩下 10, 30 Hz

allpass filter:全部都通過

頻譜平移

Spectrum Shifting 就是將輸入訊號在頻率域中平移,產生輸出訊號,也稱為頻率平移 (frequency shifting)

如果是聲音訊號經過頻譜平移,可改變聲音的頻率,因此,常被稱為音頻改變,或音高改變(Pitch Change)

回顧傅立葉轉換的第二平移定理(或頻率平移定理)

\(F\{f(t) \cdot e^{jω_0t}\} = F(ω-ω_0)\),其中 \(ω_0 > 0 \) 是平移的角頻率。 因此函數 f 的傅立葉轉換,在頻率域的平移,結果是原函數另外乘上一個複數指標函數。

本章採用的方法是在頻率域中進行平移


ex: 輸入弦波訊號的定義:\(x(t)=cos(2 \pi \cdot 10 \cdot t)\) ,且取樣頻率 \(f_s = 500Hz\),套用頻譜平移 (Spectrum Shifting) 技術,平移 20Hz 的頻率

import numpy as np
from numpy.fft import fft, fftshift, ifft, fftfreq
import matplotlib.pyplot as plt

def spectrum_shifting( x, shift, fs ):
    X = fft( x )
    N = fs
    N_half = int( fs / 2 )
    Y = np.zeros( N, dtype = 'complex' )
    for i in range( N_half ):
        if i + shift >= 0 and i + shift <= N_half:
            Y[i + shift] = X[i]
    for i in range( N_half + 1, fs ):
        if i - shift >= N_half + 1 and i - shift < N:
            Y[i - shift] = X[i]
    y = ifft( Y )
    y = y.real
    return y

def main( ):
    fs = 500
    # 輸入弦波訊號
    t = np.linspace( 0, 1, fs, endpoint = False )
    x = np.cos( 2 * np.pi * 50 * t )

    # 平移
    y = spectrum_shifting( x, -30, fs )

    f = fftshift( fftfreq( fs, 1 / fs ) )
    Xm = abs( fftshift( fft( x ) ) )
    Ym = abs( fftshift( fft( y ) ) )

    # 替換中文字型,處理 matplotlib 中文 label 問題
    # https://blog.csdn.net/gmr2453929471/article/details/78655834
    plt.rcParams['font.sans-serif'] = ['Heiti TC'] # 步驟一(替換sans-serif字型)
    plt.rcParams['axes.unicode_minus'] = False  # 步驟二(解決座標軸負數的負號顯示問題)

    plt.subplot(221)
    plt.plot( x )
    plt.xlabel( 't (second) (輸入訊號)' )
    plt.ylabel( 'Amplitude' )

    plt.subplot(222)
    plt.plot( f, Xm )
    plt.xlabel( 'f (輸入訊號頻譜)' )
    plt.ylabel( 'Magnitude' )

    plt.subplot(223)
    plt.plot( y )
    plt.xlabel( 't (second) (平移後輸出訊號)' )
    plt.ylabel( 'Amplitude' )

    plt.subplot(224)
    plt.plot( f, Ym )
    plt.xlabel( 'f (平移後輸出訊號頻譜)' )
    plt.ylabel( 'Magnitude' )

    plt.show( )

main( )

輸入訊號的頻率分量為 10Hz,平移 20 Hz 後,輸出為 30 Hz

語音的頻率域 DSP

上面的範例輸入訊號時間長度 1s,取樣率 500Hz,FFT 運算還可在有限時間內處理完畢。

但 DSP 實務上,輸入訊號的樣本數很多,例如一般來說歌曲的取樣率 44.1kHz,時間長度 3mins,單通道語音的樣本數為 \(44100*60*3=7938000\)

如果對整個輸入訊號進行 FFT,運算量過大。

因此數位語音 DSP 實務中,通常是將原始的輸入訊號,切割為多個獨立的 segments,每個 segments 包含 N 個樣本數,對每個 segment 分別進行頻率域 DSP,再重組成輸出訊號的結果。

wav ideal filtering

import numpy as np
import wave
from scipy.io.wavfile import read, write
import struct
from numpy.fft import fft, fftshift, ifft

def ideal_lowpass_filtering( x, cutoff, fs ):
    X = fft( x )
    H = np.zeros( fs )
    for i in range( -cutoff, cutoff + 1 ):
        H[i] = 1
    Y = H * X
    y = ifft( Y )
    y = y.real
    return y

def ideal_highpass_filtering( x, cutoff, fs ):
    X = fft( x )
    H = np.zeros( fs )
    for i in range( -cutoff, cutoff + 1 ):
        H[i] = 1
    H = 1 - H
    Y = H * X
    y = ifft( Y )
    y = y.real
    return y

def ideal_bandpass_filtering( x, f1, f2, fs ):
    X = fft( x )
    H = np.zeros( fs )
    for i in range( f1, f2 + 1 ):
        H[i] = 1
    for i in range( -f1, -f2 - 1, -1 ):
        H[i] = 1
    Y = H * X
    y = ifft( Y )
    y = y.real
    return y

def ideal_bandstop_filtering( x, f1, f2, fs ):
    X = fft( x )
    H = np.zeros( fs )
    for i in range( f1, f2 + 1 ):
        H[i] = 1
    for i in range( -f1, -f2 - 1, -1 ):
        H[i] = 1
    H = 1 - H
    Y = H * X
    y = ifft( Y )
    y = y.real
    return y

def ideal_allpass_filtering( x ):
    X = fft( x )
    Y = X
    y = ifft( Y )
    y = y.real
    return y

def dsp(choise, x, fs):
    y = np.zeros( x.size )
    n = int( x.size / fs ) + 1
    N = fs
    for iter in range( n ):
        xx = np.zeros( N )
        yy = np.zeros( N )
        for i in range( iter * N, ( iter + 1 ) * N ):
            if i < x.size:
                xx[i - iter * N] = x[i]

        # yy = ideal_lowpass_filtering( xx, 2000, fs )
        if choise == 0:
            yy = ideal_lowpass_filtering( xx, 2000, fs )
        elif choise == 1:
            yy = ideal_highpass_filtering( xx, 2000, fs )
        elif choise == 2:
            yy = ideal_bandpass_filtering(xx, 2000, 3000, fs)
        elif choise == 3:
            yy = ideal_bandstop_filtering(xx, 2000, 3000, fs)
        else:
            yy = ideal_allpass_filtering(xx)

        for i in range( iter * N, ( iter + 1 ) * N ):
            if i < x.size:
                y[i] = yy[i - iter * N]

    return y

def write_wav_file(y, filename, num_channels, sampwidth, fs, num_frames, comptype, compname):
    wav_file = wave.open( filename, 'w' )
    wav_file.setparams(( num_channels, sampwidth, fs, num_frames, comptype, compname ))

    for s in y:
        wav_file.writeframes( struct.pack( 'h', int ( s ) ) )

    wav_file.close( )

def main( ):
    # infile  = input( "Input File: " )
    # outfile = input( "Output File: " )
    infile = "r2d2.wav"

    # ----------------------------------------------------
    #  輸入模組
    # ----------------------------------------------------
    wav = wave.open( infile, 'rb' )
    num_channels = wav.getnchannels( )  # 通道數
    sampwidth    = wav.getsampwidth( )  # 樣本寬度
    fs           = wav.getframerate( )  # 取樣頻率(Hz)
    num_frames   = wav.getnframes( )    # 音框數 = 樣本數
    comptype     = wav.getcomptype( )   # 壓縮型態
    compname     = wav.getcompname( )   # 無壓縮
    wav.close( )

    sampling_rate, x = read( infile )   # 輸入訊號

    # ----------------------------------------------------
    #  DSP 模組
    # ----------------------------------------------------
    # y = np.zeros( x.size )
    # n = int( x.size / fs ) + 1
    # N = fs
    # for iter in range( n ):
    #   xx = np.zeros( N )
    #   yy = np.zeros( N )
    #   for i in range( iter * N, ( iter + 1 ) * N ):
    #       if i < x.size:
    #           xx[i - iter * N] = x[i]

    #   yy = ideal_lowpass_filtering( xx, 2000, fs )

    #   for i in range( iter * N, ( iter + 1 ) * N ):
    #       if i < x.size:
    #           y[i] = yy[i - iter * N]

    # ----------------------------------------------------
    #  輸出模組
    # ----------------------------------------------------
    # wav_file = wave.open( outfile, 'w' )
    # wav_file.setparams(( num_channels, sampwidth, fs, num_frames, comptype, compname ))

    # for s in y:
    #   wav_file.writeframes( struct.pack( 'h', int ( s ) ) )

    # wav_file.close( )
    y = dsp(0, x, fs)
    write_wav_file(y, "01_ideal_lowpass_filtering.wav", num_channels, sampwidth, fs, num_frames, comptype, compname)

    y = dsp(1, x, fs)
    write_wav_file(y, "02_ideal_highpass_filtering.wav", num_channels, sampwidth, fs, num_frames, comptype, compname)

    y = dsp(2, x, fs)
    write_wav_file(y, "03_ideal_bandpass_filtering.wav", num_channels, sampwidth, fs, num_frames, comptype, compname)

    y = dsp(3, x, fs)
    write_wav_file(y, "04_ideal_bandstop_filtering.wav", num_channels, sampwidth, fs, num_frames, comptype, compname)

    y = dsp(4, x, fs)
    write_wav_file(y, "05_ideal_allpass_filtering.wav", num_channels, sampwidth, fs, num_frames, comptype, compname)

main( )

頻譜平移

原始的輸入訊號,必須先切割成多個獨立的 segments,每個片段包含 N 個樣本,再分別對每個片段進行頻譜平移

import numpy as np
import wave
from scipy.io.wavfile import read, write
import struct
from numpy.fft import fft, fftshift, ifft

def spectrum_shifting( x, shift, fs ):
    X = fft( x )
    N = fs
    N_half = int( fs / 2 )
    Y = np.zeros( N, dtype = 'complex' )
    for i in range( N_half ):
        if i + shift >= 0 and i + shift <= N_half:
            Y[i + shift] = X[i]
    for i in range( N_half + 1, fs ):
        if i - shift >= N_half + 1 and i - shift < N:
            Y[i - shift] = X[i]
    y = ifft( Y )
    y = y.real
    return y

def main( ):
    # infile  = input( "Input File: " )
    # outfile = input( "Output File: " )
    infile  = "r2d2.wav"
    outfile = "r2d2_spectrum_shifting.wav"

    # ----------------------------------------------------
    #  輸入模組
    # ----------------------------------------------------
    wav = wave.open( infile, 'rb' )
    num_channels = wav.getnchannels( )  # 通道數
    sampwidth    = wav.getsampwidth( )  # 樣本寬度
    fs           = wav.getframerate( )  # 取樣頻率(Hz)
    num_frames   = wav.getnframes( )    # 音框數 = 樣本數
    comptype     = wav.getcomptype( )   # 壓縮型態
    compname     = wav.getcompname( )   # 無壓縮
    wav.close( )

    sampling_rate, x = read( infile )   # 輸入訊號

    # ----------------------------------------------------
    #  DSP 模組
    # ----------------------------------------------------
    y = np.zeros( x.size )
    n = int( x.size / fs ) + 1
    N = fs
    for iter in range( n ):
        xx = np.zeros( N )
        yy = np.zeros( N )
        for i in range( iter * N, ( iter + 1 ) * N ):
            if i < x.size:
                xx[i - iter * N] = x[i]

        yy = spectrum_shifting( xx, 500, fs )

        for i in range( iter * N, ( iter + 1 ) * N ):
            if i < x.size:
                y[i] = yy[i - iter * N]

    # ----------------------------------------------------
    #  輸出模組
    # ----------------------------------------------------
    wav_file = wave.open( outfile, 'w' )
    wav_file.setparams(( num_channels, sampwidth, fs, num_frames, comptype, compname ))

    for s in y:
        wav_file.writeframes( struct.pack( 'h', int ( s ) ) )

    wav_file.close( )

main( )

References

數位訊號處理:Python程式實作(附範例光碟)(第二版)