2020/11/30

數位訊號雜訊 in python

雜訊生成

\(y[n]=x[n]+η[n]\) 其中 η[n] 稱為雜訊 Noise

DSP 中,如果資料分布情形或統計參數(例如標準差、平均值),不會隨著時間改變,則該訊號為平穩訊號 Stationary Signals,反之稱為非平穩訊號。例如馬達產生的聲音是平穩訊號,一般語音是非平穩訊號

均勻雜訊 Uniform Noise

以 uniform distribution function 產生隨機亂數資料,因為是均勻分佈,平均值 𝜇 為 0,雜訊連續樣本之間互相獨立,沒有相關。在 DSP 這種雜訊稱為 White Noise

均勻分佈函數:

給定隨機變數 random variable z,Uniform Distribution Function 定義為

\(p(z) = \left\{ \begin{array}{ll} 1/2 & \mbox{if -1≤z≤1} \\ 0 & \mbox{otherwise} \end{array} \right.\) ,其中 z 介於 -1~1 之間

均勻分佈的機率密度函數,必須滿足總機率值 (積分) 為 1。曲線下的面積為 1

\(\int_{-∞}^{∞} p(z) {\rm d}z = 1\)

高斯雜訊 Gaussian Noise

高斯分佈函數定義

\(p(z) = \frac{1}{\sqrt{2 \pi σ^2}} e^{(z-μ)^2/(2σ^2)}\)

其中 μ 為平均,σ 是標準差

噪音的隨機亂數資料以高斯分佈函數也就是常態分佈函數 Normal Distribution Function 產生,其平均值 𝜇 為 0,標準差 𝜎 為 1,且滿足 \(\int_{-∞}^{∞}p(z)dz = 1\),就是函數曲線下的面積總和為 1(所有機率的總和為 1)

雜訊連續樣本之間互相獨立,沒有自相關 autocorrelation。在 DSP 這種雜訊稱為 Gaussian White Noise

布朗尼雜訊 Brownian Noise

根據布朗尼運動 Brownian Motion 產生的亂數資料,也稱為 Random Walk Noise

Brownian Noise 也稱為 Brown Noise,因為 Brownian Motion 的發現者 Robert Brown 而這樣命名。

Brownian Noise 的每個 sample,取自先前所有 samples 的總和(積分)。

程式中產生 Brownian Noise 後,會再針對數值做正規化 Normalization,讓結果介於 -1~1 之間

Brownian Noise 的變動幅度比均勻雜訊緩慢,連續的 samples 之間有相關性,結果會造成數位訊號的振幅微幅改變的狀況,也就是發生訊號飄移的狀況

脈衝雜訊 Impulse Noise

瞬間的脈衝產生的雜訊,當訊號受到電磁干擾、或是唱片(CD)有刮痕,會產生 Impulse Noise。可用振幅 Amplitude 與發生機率 Probability 來模擬產生 Impulse Noise

以 python 程式碼產生 noise

需要的套件

pip install numpy
pip install scipy

程式裡面有以下這些雜訊的產生方式

import numpy as np
import scipy.signal as signal
import wave
import struct
import sys

class WavFileParams:
    def __init__(self, filename, amplitude, frequency, duration, fs, num_channels, sampwidth, comptype, compname):
        self.filename = filename
        self.amplitude = amplitude                      # 振幅
        self.frequency = frequency                      # 頻率(Hz)
        self.duration = duration                        # 時間長度(秒)
        self.fs = fs                                    # 取樣頻率(Hz)
        self.num_samples = self.duration * self.fs      # 樣本數

        self.num_channels = num_channels                # 通道數
        self.sampwidth = sampwidth                      # 樣本寬度
        self.num_frames = self.num_samples              # 音框數 = 樣本數
        self.comptype = comptype                        # 壓縮型態
        self.compname = compname                        # 壓縮名稱

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

    for s in digital_signel_points :
        # 以 struct.pack 將 int 資料轉換為 16 bits
        wav_file.writeframes( struct.pack( 'h', int ( s ) ) )

    wav_file.close( )

def uniform_noise():
    wav_params = WavFileParams(
        filename = "uniform_noise.wav",     # 檔案名稱

        amplitude = 30000,          # 振幅
        frequency = 100,            # 頻率(Hz)
        duration = 1,               # 時間長度(秒)
        fs = 44100,                 # 取樣頻率(Hz)
        # num_samples = duration * fs,  # 樣本數

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

    t = np.linspace( 0, wav_params.duration, wav_params.num_samples, endpoint = False );
    # 原始訊號
    x = wav_params.amplitude * np.cos( 2 * np.pi * wav_params.frequency * t );
    # uniform noise, 在 -1000 ~ 1000 之間
    noise = np.random.uniform( -1000, 1000, wav_params.num_samples )
    y = x + noise

    write_wav_file(wav_params, y);

def gaussian_noise():
    wav_params = WavFileParams(
        filename = "gaussian_noise.wav",        # 檔案名稱

        amplitude = 30000,          # 振幅
        frequency = 100,            # 頻率(Hz)
        duration = 1,               # 時間長度(秒)
        fs = 44100,                 # 取樣頻率(Hz)
        # num_samples = duration * fs,  # 樣本數

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

    t = np.linspace( 0, wav_params.duration, wav_params.num_samples, endpoint = False );
    # 原始訊號
    x = wav_params.amplitude * np.cos( 2 * np.pi * wav_params.frequency * t );
    # 雜訊(高斯分佈), 在 0 ~ 1000 之間
    noise = np.random.normal( 0, 1000, wav_params.num_samples )
    y = x + noise

    # 為避免相加後超過 16 bits 資料長度的值,選擇使用比較小的振幅,且在相加後,用 np.clip 抑制數值範圍
    y_clip = np.clip(y, -32768, 32767)

    write_wav_file(wav_params, y_clip);

def brownian_noise():
    wav_params = WavFileParams(
        filename = "brownian_noise.wav",        # 檔案名稱

        amplitude = 30000,          # 振幅
        frequency = 100,            # 頻率(Hz)
        duration = 1,               # 時間長度(秒)
        fs = 44100,                 # 取樣頻率(Hz)
        # num_samples = duration * fs,  # 樣本數

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

    t = np.linspace( 0, wav_params.duration, wav_params.num_samples, endpoint = False );
    # 原始訊號
    x = wav_params.amplitude * np.cos( 2 * np.pi * wav_params.frequency * t );

    # 布朗尼雜訊 brownian_noise, 在 -1 ~ 1 之間
    n1 = np.random.uniform( -1, 1, wav_params.num_samples )
    # 累加
    ns = np.cumsum( n1 )
    # 正規化
    mean = np.mean( ns )
    max = np.max( np.absolute( ns - mean ) )
    noise = (( ns - mean ) / max)
    # 把最後正規化的數值,放大 10000 倍,強化 noise
    noise = noise * 10000.0

    y = x + noise

    # 為避免相加後超過 16 bits 資料長度的值,選擇使用比較小的振幅,且在相加後,用 np.clip 抑制數值範圍
    y_clip = np.clip(y, -32768, 32767)

    write_wav_file(wav_params, y_clip);

def impulse_noise():
    wav_params = WavFileParams(
        filename = "impulse_noise.wav",     # 檔案名稱

        amplitude = 30000,          # 振幅
        frequency = 100,            # 頻率(Hz)
        duration = 1,               # 時間長度(秒)
        fs = 44100,                 # 取樣頻率(Hz)
        # num_samples = duration * fs,  # 樣本數

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

    # impulse noise 發生機率 (%)
    probability = 5

    t = np.linspace( 0, wav_params.duration, wav_params.num_samples, endpoint = False );
    # 原始訊號
    x = wav_params.amplitude * np.cos( 2 * np.pi * wav_params.frequency * t );

    noise = np.zeros( x.size )                          # 脈衝雜訊
    for i in range( x.size ):
        p1 = np.random.uniform( 0, 1 )
        if p1 < probability / 100:
            p2 = np.random.uniform( 0, 1 )
            if p2 < 0.5:
                noise[i] = wav_params.amplitude
            else:
                noise[i] = -wav_params.amplitude

    y = x + noise

    # 為避免相加後超過 16 bits 資料長度的值,選擇使用比較小的振幅,且在相加後,用 np.clip 抑制數值範圍
    y_clip = np.clip(y, -32768, 32767)

    write_wav_file(wav_params, y_clip);

def main():
    uniform_noise()
    gaussian_noise()
    brownian_noise()
    impulse_noise()

if __name__ == "__main__":
    sys.exit(main())

以下分別是產生出來的 wav 的波形,第三個 Brownian Noise 可看出有上下飄移的狀況

訊號雜訊比

Signal-to-Noise Ratio 也就是 SNR 定義為

\( SNR = \frac{P_{signal}}{P_{noise}} \) ,\(P_{signal}\) 是訊號功率,\(P_{noise}\) 是雜訊功率

SNR 也可以定義為

\( SNR = (\frac{A_{signal}}{A_{noise}})^2 \) ,\(A_{signal}\) 是訊號振幅,\(A_{noise}\) 是雜訊振幅

SNR 通常用分貝 Decibels dB 為單位,故定義為

\( SNR = 10log_{10}(\frac{P_{signal}}{P_{noise}}) dB\) 或

\( SNR = 10log_{10}(\frac{A_{signal}}{A_{noise}})^2 dB = 20log_{10}(\frac{A_{signal}}{A_{noise}}) dB \)

如果 SNR 值越高,表示訊號的品質越好, SNR 越低則訊號受到雜訊干擾越嚴重

References

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

沒有留言:

張貼留言