2021/04/12

mongodb 帳號驗證權限

資料庫的權限有四種

  1. readAnyDatabase 任何資料庫的唯讀權限
  2. userAdminAnyDatabase 任何資料庫的讀寫權限
  3. userAdminAnyDatabase 任何資料庫用戶的管理權限
  4. dbAdminAnyDatabase 任何資料庫的管理權限

以往在啟動 mongodb 時,就直接指定 mongod.conf

mongod -f $MONGODB_HOME/conf/mongod.conf

這時候可用 mongo client 連接資料庫

./mongo

查看有沒有 users

db.system.users.find()

查看 user

show users

建立 admin 資料庫的管理者帳號

use admin

db.createUser({
  user : "root",
  pwd : "password",
  roles : [
    "clusterAdmin",
    "dbAdminAnyDatabase",
    "userAdminAnyDatabase",
    "readWriteAnyDatabase"
  ]
})

重新啟動 mognodb (加上 --auth 參數)

mongod --auth -f $MONGODB_HOME/conf/mongod.conf

再次查看 users

> use admin
> db.system.users.find()
Error: error: {
    "ok" : 0,
    "errmsg" : "command find requires authentication",
    "code" : 13,
    "codeName" : "Unauthorized"
}
> show users
2020-09-02T11:36:48.994+0800 E QUERY    [js] Error: command usersInfo requires authentication :
_getErrorWithCode@src/mongo/shell/utils.js:25:13
DB.prototype.getUsers@src/mongo/shell/db.js:1763:1
shellHelper.show@src/mongo/shell/utils.js:859:9
shellHelper@src/mongo/shell/utils.js:766:15
@(shellhelp2):1:1

驗證後,就可以查看 users

> use admin
switched to db admin
> db.auth("root", "password")
1
> show users
{
    "_id" : "admin.root",
    "user" : "root",
    "db" : "admin",
    "roles" : [
        {
            "role" : "clusterAdmin",
            "db" : "admin"
        },
        {
            "role" : "dbAdminAnyDatabase",
            "db" : "admin"
        },
        {
            "role" : "userAdminAnyDatabase",
            "db" : "admin"
        },
        {
            "role" : "readWriteAnyDatabase",
            "db" : "admin"
        }
    ],
    "mechanisms" : [
        "SCRAM-SHA-1",
        "SCRAM-SHA-256"
    ]
}

針對資料庫,建立管理者

use admin

db.createUser({
  user : "maxkit",
  pwd : "password",
  roles : [
    {role:"readWrite", db:"larzio"}
  ]
})

針對資料庫,建立一般使用者

use larzio

db.createUser({
  user : "maxkit",
  pwd : "max168kit",
  roles : [
    {role:"readWrite", db:"larzio"}
  ]
})

刪除帳號

use admin

db.dropUser("maxkit")

關閉 mognodb

mongo 127.0.0.1:27017 -u root -p 'password' --authenticationDatabase 'admin' --eval "db.getSiblingDB('admin').shutdownServer()"

References

Mongodb 創建管理員帳號與普通帳號

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程式實作(附範例光碟)(第二版)

2021/02/22

頻率響應

DSP (ex: 音響、通訊設備) 的目的都是希望達到無失真的傳遞,在輸出端重現原始訊號,因此,理想的訊號處理系統,在不同頻率範圍,需表現均勻的強度響應。頻率響應 Frequency Response 是用來衡量系統在特定頻率範圍的操作特性。

音響設備可用頻率響應來衡量品質,ex: 理想的喇叭,在人類聽力範圍中 20Hz~20kHz,應該具備等值的頻率響應,藉以達到原音重現。

頻率響應的定義:

系統對於輸入訊號在特定頻率範圍的響應情形,可定義為頻率的函式,響應則包含強度 Magnitude 與相位角 Phase Angle 等等量化數據,通常是以頻譜的方式呈現,用來表示系統的操作特性。


典型的 DSP (LTI) 系統, x[n] 為輸入訊號,經過 h[n] 脈衝響應(impulse response, 也稱為濾波器 filter),得到 y[n] 輸出訊號

如果 DSP 牽涉卷積運算 \(y[n]=h[n]*x[n]\)

根據卷積定理 Covolution Theorem 則

\(F\{y[n]\} = F\{h[n]*x[n]\} = F\{h[n]\} \cdot F\{x[n]\}\)

或是以離散時間傅立葉轉換可表示成:

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

因此系統的頻率響應為:

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

根據 z 轉換,系統的頻率響應跟轉換函式的關係,可表示為:

\(H(e^{jω}) = H(z)|_{z=e^{jω}}\)

因此系統的頻率響應,可表示為

\(H(z) = \frac{Y(z)}{X(z)}\)

其中 X(z), Y(z) 分別為輸入訊號與輸出訊號的 z 轉換結果


DSP 系統的頻率響應 定義為:

\(H(e^{jω}) = \sum_{n=-∞}^{∞} h[n]e^{-jωn}\)

在 LTI 系統,x[n] 是輸入訊號,y[n] 是輸出訊號,卷積運算為:

\(y[n]=\sum_{k=-∞}^{∞}h[k] \cdot x[n-k]\)

假設輸入訊號為複數指數訊號 \(x[n]=e^{jωn}\)

則輸出訊號為 \(y[n]=\sum_{k=-∞}^{∞}h[k] \cdot e^{jω(n-k)} = (\sum_{k=-∞}^{∞}h[k] \cdot e^{jωk})e^{jωn}\)

可改寫為 \(y[n]=H(e^{jω})e^{jωn}\)

因此,可得到以下公式:

\(H(e^{jω}) = \sum_{n=-∞}^{∞}h[n]e^{-jωn}\)

這個公式稱為 DSP 的頻率響應 frequency response

頻率響應就是對脈衝響應 h[n] 求離散時間傅立葉轉換 DTFT 的結果

因為 DTFT 的結果通常是複數,因此經常取其強度 magnitude,與相位角 phase angle

\(|H(e^{jω})|\) 與 \(ang\{H(e^{jω})\}\)

結果都是實數。以圖形表示就是強度頻譜 magnitude spectrum 與相位頻譜 phase spectrum,這是在頻率域的操作特性

頻率響應也常以分貝 decibels, dB 為單位表示,稱為系統的增益 Gain

\(G(ω) = 20 \cdot log_{10}|H(e^{jω}| (dB)\)

濾波器分類

根據系統(濾波器)的頻率響應,可分成以下幾種:

  • 低通濾波器 Lowpass Filter

    使得低頻範圍的訊號通過,抑制高頻範圍的訊號,頻率的閥值稱為截止頻率(Cutoff Frequency) ,定義為 \(f_c\),對應的截止角頻率 \(ω_c\)

  • 高通濾波器 Highpass Filter

    低通濾波器的相反,讓高頻範圍訊號通過,抑制低頻訊號

  • 帶通濾波器 Bandpass Filter

    讓頻率介於 \(f_1 ~ f_2\) 或角頻率 \(ω_1~ω_2\) 之間的訊號通過,抑制其他頻率範圍的訊號

  • 帶阻濾波器 Bandstop Filter

    帶通濾波器的相反,抑制頻率範圍落在 \(f_1 ~ f_2\) 或角頻率 \(ω_1~ω_2\) 之間的訊號,使得其他頻率範圍的訊號通過

  • 陷波濾波器 Notch Filter

    是帶阻濾波器的一種,僅抑制某個特定頻率,ex: 60Hz 的干擾雜訊。陷波濾波器主要是使得某特定頻率的訊號變為 0,因此也稱為歸零濾波器 Nulling Filter

  • 梳狀濾波器 Comb Filter

    Comb Filter 其頻率響應的形狀跟梳子一樣,主要是同時抑制多個頻率範圍,頻率與頻率之間,通常是設定為整數倍數

  • 全通濾波器 All-Pass Filter

    讓所有頻率範圍都通過,通常會調整訊號在某些頻率範圍的相位移

濾波器也可根據訊號的型態分成 類比濾波器 Analog Filter,與數位濾波器 Digital Filter

頻率響應範例

介紹幾種經典的頻率響應

理想濾波器

濾波器的用途是讓某些特定範圍的頻率分量通過,抑制其他範圍的頻率分量,為了不造成訊號失真,頻率通過的頻率分量的頻率響應值是設定為 1,抑制的頻率分量,頻率響應值設定為 0。

使得頻率分量通過的頻率範圍稱為通過頻帶 (Passband),抑制的頻率範圍稱為截止頻帶 (Stopband)

  • 理想低通濾波器 Ideal Lowpass Filter

    轉換函式:\( H_{LP}(e^{jω}) = \left\{ \begin{array}{ll} 1 & \mbox{if 0 ≤ ω ≤ \)ω_c\(} \\ 0 & \mbox{otherwise} \end{array} \right.\)

    在此只定義單邊 Single-Sided 的頻率響應,ω 介於 0 ~ π 之間,且 \(ω_c=2 \pi f_c\),其中 \(f_c\) 稱為截止頻率 Cutoff Frequency

  • 理想高通濾波器 Ideal Highpass Filter

    轉換函式:\( H_{HP}(e^{jω}) = \left\{ \begin{array}{ll} 1 & \mbox{if \)ω_c\( < ω ≤ π} \\ 0 & \mbox{otherwise} \end{array} \right.\)

    低通與高通濾波器的關係為: \(H_{HP}(e^{jω})=1- H_{LP}(e^{jω})\)

  • 理想帶通濾波器 Ideal Bandpass Filter

    轉換函式:\( H_{BP}(e^{jω}) = \left\{ \begin{array}{ll} 1 & \mbox{if \)ω1\( ≤ ω ≤ \)ω2\(} \\ 0 & \mbox{otherwise} \end{array} \right.\)

    允許頻率範圍 \(ω_1 ≤ ω ≤ ω_2\) 或頻率 \(f_1 ≤ f ≤ f_2\) 的訊號通過,\(f_1, f_2\) 稱為截止頻率

  • 理想帶阻濾波器 Ideal Bandstop Filter

    轉換函式:\( H_{BS}(e^{jω}) = \left\{ \begin{array}{ll} 1 & \mbox{if ω < \)ω1\( or ω > \)ω2\(} \\ 0 & \mbox{otherwise} \end{array} \right.\)

    帶通與帶阻濾波器的關係為: \(H_{BS}(e^{jω})=1- H_{BP}(e^{jω})\)


ex: Ideal Lowpass Filter定義如下,求其在離散時間域的脈衝響應 Impulse Response

​ \( H_{LP}(e^{jω}) = \left\{ \begin{array}{ll} ​ 1 & \mbox{if 0 ≤ ω ≤ \)ω_c\(} \\ ​ 0 & \mbox{otherwise} ​ \end{array} \right.\)

根據 Inverse DTFT 公式

\(x[n] = \frac{1}{2π} \int_{-∞}^{∞}X(e^{jωn})e^{jωn}dω\)

Ideal Lowpass Filter 在離散時間域的脈衝響應為:

\(h_{LP}[n] = \frac{1}{2π} \int_{-∞}^{∞}H_{LP}(e^{jωn})e^{jωn}dω \\ = \frac{1}{2π} \int_{-ω_c}^{ω_c}e^{jωn}dω \\ = \frac{1}{2π} [\frac{1}{jn}e^{jωn}]_{-ω_c}^{ω_c} \\ = \frac{1}{2π} [\frac{e^{jω_cn}}{jn}- \frac{e^{-jω_cn}}{jn} ] \\ = \frac{sin(ω_cn)}{πn} \quad\quad \mbox{-∞<n<∞}\)

當 n=0,使用羅必達規則 L'Hopital's Rule

\(h_{LP}[0] = \lim_{n→0} \frac{sin(ω_cn)}{πn} \\ = \lim_{n→0} \frac{ \frac{d}{dn} sin(ω_cn)}{ \frac{d}{dn} (πn) } \\ = \lim_{n→0} \frac{ cos(ω_cn) \cdot ω_c}{ π } \\ = \frac{ ω_c }{ π }\)

因為 \(-∞ < n < ∞\),無法在離散時間域中,定義有限長度的脈衝響應,藉以實現頻率域中的 ideal lowpass filter。脈衝響應在計算時有用到過去、現在、未來的數位訊號,不是因果系統。脈衝響應不符合絕對可加總 Absolutely Summable 原則,系統也不具 BIBO 穩定性。

用 python 處理離散時間域的脈衝響應時,因為離散序列是有限長度,以近似方式接近 ideal lowpass filter。原則上採用的濾波器大小,是使用奇數

當 filter 大小為 5時,-2<n<2

當 filter 大小為 11 時,-5<n<5

ideal lowpass filter 的脈衝響應,具有 sinc 函數的特性,向左右兩邊無限延伸

以下範例選擇 \(ω_c=\frac{\pi}{2}\),因此當 n=0,\(\frac{ω_c}{\pi}=0.5\) ,可調整 \(ω_c\) 觀察脈衝響應的結果

import numpy as np
import matplotlib.pyplot as plt

def plot(filter_size, plotpos):
    filter_half = int( filter_size / 2 )
    wc = np.pi / 2

    na = np.arange( -filter_half, filter_half + 1 ) # 定義 n 陣列
    h = np.zeros( filter_size )                     # 計算脈衝響應
    for n in na:
        if n == 0:
            h[n+filter_half] =  wc/np.pi
        else:
            h[n+filter_half] = np.sin( wc * n ) / ( np.pi * n )

    plt.subplot(plotpos)
    plt.stem( na, h, use_line_collection=True )
    plt.xlabel( 'n (filter_size='+ str(filter_size) +')')
    plt.ylabel( 'h[n]' )

plot(5, 221)
plot(11, 222)
plot(21, 223)
plot(31, 224)

plt.show( )

以下選擇 \(ω_c=\frac{\pi}{10}\),因此當 n=0,\(\frac{ω_c}{\pi}=1/10\)


由於有限長度的脈衝響應,無法提供 ideal filter 的操作特性,在 DSP 實際應用時,通常是讓頻率響應在通過頻帶與截止頻帶之間,由 1 逐漸降為 0,且在通過頻帶或截止頻帶內,頻率響應允許有些微變動。

剛剛是求濾波器大小 5, 11,21, 31 的脈衝響應,現在改求頻率響應。

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


def plot(filter_size, plotpos):
    filter_half = int( filter_size / 2 )
    wc = np.pi / 2

    na = np.arange( -filter_half, filter_half + 1 ) # 定義 n 陣列
    h = np.zeros( filter_size )                     # 計算脈衝響應
    for n in na:
        if n == 0:
            h[n+filter_half] =  wc/np.pi
        else:
            h[n+filter_half] = np.sin( wc * n ) / ( np.pi * n )

    w, H = signal.freqz( h )
    mag = abs( H )

    plt.subplot(plotpos)
    plt.plot( w, mag )
    plt.xlabel( r'$\omega$'+' (filter_size='+ str(filter_size) +')' )
    plt.ylabel( 'Magnitude' )


plot(5, 221)
plot(11, 222)
plot(21, 223)
plot(31, 224)

plt.show( )

結果發現,若 filter size 越大,或脈衝響應樣本越多,越接近理想的 lowpass filter。此外可注意到頻率響應在通過頻帶與截止頻帶之間,出現了 Gibbs 現象。

平均濾波器

average filter 定義為

\(h[n] = \frac{1}{M}{1,1,1,...,1}, n=0,1,...,M-1\) ,M為濾波器大小 filter size

平均濾波器的頻率響應,可根據其 DTFT 得來

\(H(e^{jω}) = \sum_{n=-∞}^{∞}h[n]e^{-jωn} \\ = \frac{1}{M} \sum_{n=0}^{M-1}e^{-jωn} \\ = \frac{1}{M} ( \sum_{n=0}^{∞}e^{-jωn} - \sum_{n=M}^{∞}e^{-jωn}) \\ = \frac{1}{M} ( \sum_{n=0}^{∞}e^{-jωn} ) (1 - e^{-jωM}) \\ = \frac{1}{M} \frac{1 - e^{-jωM}}{1 - e^{-jω}} \\ = \frac{1}{M} \frac{sin(Mω/2)}{sin(ω/2)} e^{-j(M-1)ω/2} \)

取其強度 magnitude

\(|H(e^{jω})| = |\frac{sin(Mω/2)}{sin(ω/2)}|\)

以圖形表示,就是平均濾波器的頻率響應


ex: 若平均濾波器定義 \(h[n]=\frac{1}{M}{1,1,...,1}, n=0,1,...,M-1\)

當 M=5, 15,求其頻率響應

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

def average_filter(filter_size):
    h = np.ones( filter_size ) / filter_size
    return h

def plot(filter_size, plotpos):
    h = average_filter(filter_size)

    w, H = signal.freqz( h )
    mag = abs( H )

    plt.subplot(plotpos)
    plt.plot( w, mag )
    plt.xlabel( r'$\omega$'+' (filter_size='+ str(filter_size) +')' )
    plt.ylabel( 'Magnitude' )

plot(5, 121)
plot(15, 122)

plt.show( )

平均濾波器有 lowpass filter 的特性,當 size 越小,通過頻帶越寬。但平均濾波器有允許少量高頻分量通過

高斯濾波器

Gaussian Filter 可定義為

\(g[n]=e^{-n^2/2σ^2}\) ,其中 σ 為標準差

chap9: 高斯函數的傅立葉轉換,形成另一個高斯函數。因此高斯函數的頻率響應,可用另一個高斯函數表示

ex: Gaussian Filter 定義為 \(g[n]=e^{-n^2/2σ^2}\) ,當 σ=1, 3,求其頻率響應

note: 程式中有對係數進行正規化,讓 \(\sum_{n}g[n]=1\)

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

def plot(sigma, plotpos):
    filter_size = int( 6 * sigma + 1 )              # 濾波器大小
    gauss = signal.gaussian( filter_size, sigma )   # 濾波器係數
    sum = np.sum( gauss )                           # 正規化
    gauss = gauss / sum

    w, H = signal.freqz( gauss )
    mag = abs( H )

    plt.subplot(plotpos)
    plt.plot( w, mag )
    plt.xlabel( r'$\omega$'+' (sigma='+ str(sigma) +')' )
    plt.ylabel( 'Magnitude' )

plot(sigma=1, plotpos=121)
plot(sigma=3, plotpos=122)

plt.show( )

高斯濾波器有 lowpass filter 的特性,標準差越小,通過頻帶較寬。在高頻範圍的抑制效果,比平均濾波器好。

FIR 濾波器

一階 FIR 濾波器有 lowpass, highpass filter 兩種

FIR filters 裡面最簡單的是移動平均濾波器 Moving Averge Filter,濾波器大小為 2

最簡單的 FIR Lowpass Filter 的一階轉換函式定義為:

\(H_0(z)=\frac{1}{2}(1+z^{-1})\)

最簡單的 FIR Highpass Filter 的一階轉換函式定義為:

\(H_1(z)=\frac{1}{2}(1-z^{-1})\)

因為 \(H_1(z) = 1 - H_0(z) = 1-\frac{1}{2}(1+z^{-1}) = \frac{1}{2}(1-z^{-1})\)

ex: 分別求上面兩個 FIR filters 的頻率響應

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

b = np.array( [ 0.5, 0.5 ] )    # 低通濾波器
w, H0 = signal.freqz( b )
H0 = abs( H0 )

b = np.array( [ 0.5, -0.5 ] )   # 高通濾波器
w, H1 = signal.freqz( b )
H1 = abs( H1 )

plt.subplot(121)
plt.plot( w, H0 )
plt.xlabel( r'$\omega$' +' (lowpass)' )
plt.ylabel( 'Magnitude' )

plt.subplot(122)
plt.plot( w, H1 )
plt.xlabel( r'$\omega$' +' (highpass)' )
plt.ylabel( 'Magnitude' )

plt.show( )

根據 Lowpass Filter 的頻率響應,強度的最大/最小值分別是 1 與 0。頻率響應也是從 1 逐漸降到 0。

在討論 Lowpass Filter 時,當截止頻率 \(ω=ω_0\),強度降為

\(|H_0(e^{jω_c})| = \frac{1}{\sqrt{2}}|H_0(e^{j0})|\)

如果以分貝表示則為

\(20 log_{10}|H_0(e^{jω_c})| = 20 log_{10}( \frac{1}{\sqrt{2}}|H_0(e^{j0})| ) \\ = 20 log_{10}|H_0(e^{j0})| - 20 log_{10}\sqrt{2} \\ = 0 - 3.0103 ≈ -3dB\)

\(ω_c\)或 \(f_c\) 稱為 3dB 截止頻率 (3dB Cutoff Frequency),通常被視為是通過頻帶 (Passband) 的邊緣。因此,\(ω_c\)或 \(f_c\) 也常被稱為 通過頻帶邊緣頻率(Passband Edge Frequency)


串接多個 FIR Filter,可用來近似 ideal lowpass filter

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

# order = eval( input( "Enter number of cascade FIR filters: " ) )

def plot(order, plotpos):
    b = np.array( [ 0.5, 0.5 ] )
    w, H = signal.freqz( b )

    for i in range( order - 1 ):
        H *= H

    H0 = abs( H )

    plt.subplot(plotpos)
    plt.plot( w, H0 )
    plt.xlabel( r'$\omega$' +'(order='+str(order)+')' )
    plt.ylabel( 'Magnitude' )

plot(1, 221)
plot(2, 222)
plot(3, 223)
plot(4, 224)

plt.show( )

IIR 濾波器

  • IIR Lowpass Filter

    一階 IIR Lowpass Filter 轉換函式定義為 \(H_0(z)= \frac{1-α}{2} \frac{1+z^{-1}}{1-αz^{-1}}\) ,其中 \(|α|<1\) 為穩定的 IIR Filter

  • IIR Highpass Filter

    一階 IIR Highpass Filter 轉換函式定義為 \(H_1(z)= \frac{1+α}{2} \frac{1-z^{-1}}{1-αz^{-1}}\) ,其中 \(|α|<1\) 為穩定的 IIR Filter

ex: 一階 IIR Lowpass/Highpass Filter,當 α = 0.2, 0.4, 0.6, 0.8 時,求其頻率響應

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

# alpha = eval( input( "Please enter alpha (< 1): " ) )

def iir_lowpass(alpha, plotpos):
    b = np.array( [ 1 - alpha, 1 - alpha ] )
    a = np.array( [ 2, -2 * alpha ] )
    w, H = signal.freqz( b, a )
    H0 = abs( H )

    plt.subplot(plotpos)
    plt.plot( w, H0 )
    plt.xlabel( r'$\omega$' +' (alpha='+str(alpha)+')' )
    plt.ylabel( 'Magnitude' )

def iir_highpass(alpha, plotpos):
    b = np.array( [ 1 + alpha, - ( 1 + alpha ) ] )
    a = np.array( [ 2, -2 * alpha ] )
    w, H = signal.freqz( b, a )
    H1 = abs( H )

    plt.subplot(plotpos)
    plt.plot( w, H1 )
    plt.xlabel( r'$\omega$' +' (alpha='+str(alpha)+')' )
    plt.ylabel( 'Magnitude' )

plt.figure( 1 )
iir_lowpass(0.2, 221)
iir_lowpass(0.4, 222)
iir_lowpass(0.6, 223)
iir_lowpass(0.8, 224)

plt.figure( 2 )
iir_highpass(0.2, 221)
iir_highpass(0.4, 222)
iir_highpass(0.6, 223)
iir_highpass(0.8, 224)

plt.show( )

IIR Lowpass Filters

IIR Highpass Filters

  • IIR 帶通濾波器

    二階 IIR Bandpass Filter 的轉換函式定義為

    \(H_{BP}(z) = \frac{1-α}{2} \frac{1-z^{-2}}{1-β(1+α)z^{-1}+αz^{-2}}\) ,其中 α, β 都小於 1

ex1: 當 α = 0.2, 0.5, 0.8 且 β=0.5 時,求其頻率響應

ex2: 當 α = 0.5 且 β=0.2, 0.5, 0.8 時,求其頻率響應

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

def IIR_bandpass_frequency_response(alpha, beta, plotpos):
    b = np.array( [ 1 - alpha, 0, -( 1 - alpha ) ] )
    a = np.array( [ 2, -2 * beta * ( 1 + alpha ), 2 * alpha ] )
    w, H = signal.freqz( b, a )
    H1 = abs( H )

    plt.subplot(plotpos)

    plt.plot( w, H1, '-', label = r'$\alpha$ = '+str(alpha)+','+r'$\beta$ = '+str(beta) )
    plt.legend( loc = 'upper right' )
    plt.xlabel( r'$\omega$' )
    plt.ylabel( 'Magnitude' )


IIR_bandpass_frequency_response(0.2, 0.5, 121)
IIR_bandpass_frequency_response(0.5, 0.5, 121)
IIR_bandpass_frequency_response(0.8, 0.5, 121)

IIR_bandpass_frequency_response(0.5, 0.2, 122)
IIR_bandpass_frequency_response(0.5, 0.5, 122)
IIR_bandpass_frequency_response(0.5, 0.8, 122)

plt.show( )

α 參數用來調整頻寬,β 參數用來調整頻率中心

梳狀濾波器

  • Comb Lowpass Filter

    簡易梳狀低通濾波器的轉換函式,可定義為

    \(H_{CLP}(z) = \frac{1}{2}(1+z^{-M})\),M 為正整數

  • Comb Highpass Filter

    簡易梳狀高通濾波器的轉換函式,可定義為

    \(H_{CHP}(z) = \frac{1}{2}(1-z^{-M})\) ,M 為正整數

ex: 當 M=10,求頻率響應

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

# M = eval( input( "Please enter M for the Comb Filter: " ) )

def comb_lowpass_filter_frequency_response(M, plotpos):

    b = np.zeros( M + 1 )           # 梳狀低通濾波器
    b[0] = 0.5
    b[M] = 0.5
    w, H = signal.freqz( b, 1 )
    H0 = abs( H )

    plt.subplot(plotpos)
    plt.plot( w, H0 )
    plt.xlabel( r'$\omega$'+' (comb lowpass)' )
    plt.ylabel( 'Magnitude' )

def comb_highpass_filter_frequency_response(M, plotpos):

    b = np.zeros( M + 1 )
    b[0] = 0.5
    b[M] = -0.5                     # 梳狀高通濾波器
    w, H = signal.freqz( b, 1 )
    H1 = abs( H )

    plt.subplot(plotpos)
    plt.plot( w, H1 )
    plt.xlabel( r'$\omega$'+' (comb highpass)' )
    plt.ylabel( 'Magnitude' )

M = 10
comb_lowpass_filter_frequency_response(M, 121)
comb_highpass_filter_frequency_response(M, 122)

plt.show( )

References

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