數位訊號生成 in python

週期性訊號

\(x(t)=x(t+T)\) ，其中 T 為週期，是固定值

弦波 Sinusoids

\(x[n] = x(nT_s) = A cos(2 \pi f nT_s)= A cos(2 \pi fn/f_s)\)

ex: \(x(t)=cos(2 \pi \cdot 5 t)\) 其中 A=1, f=5Hz, 時間 0~1s

import numpy as np
import matplotlib.pyplot as plt

t = np.linspace( 0, 1, 1000, endpoint = False ) # 定義時間陣列
x = np.cos( 2 * np.pi * 5 * t )                 # 產生弦波

plt.plot( t, x )                                # 繪圖
plt.xlabel( 't (second)' )
plt.ylabel( 'Amplitude' )
plt.axis( [ 0, 1, -1.2, 1.2 ] )

plt.show( )

方波 Square

\(x(t) = A \cdot sgn(sin(ωt))\) 或 \(x(t)=A \cdot sgn(sin(2 \pi ft))\)

其中 A 為振幅，ω 是角頻率 Angular Frequency，f 為頻率 frequency

sgn 是符號函數 \(sgn(x) = \left\{ \begin{array}{ll} 1 & \mbox{if x>0} \\ 0 & \mbox{if x = 0} \\ -1 & \mbox{if x<0} \end{array} \right.\)

可用 SciPy 的 signal.square 產生方波

ex: A =1, f=5Hz, t=0~1 的方波

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

t = np.linspace( 0, 1, 1000 )           # 定義時間陣列
x = signal.square( 2 * np.pi * 5 * t )  # 產生方波

plt.plot( t, x )                        # 繪圖
plt.xlabel( 't (second)' )
plt.ylabel( 'Amplitude' )
plt.axis( [ 0, 1, -1.2, 1.2 ] )

plt.show( )

鋸齒波 Sawtooh Wave

A = 1, f=5Hz, t=0~1s 的鋸齒波

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

t = np.linspace( 0, 1, 1000 )            # 定義時間陣列
x = signal.sawtooth( 2 * np.pi * 5 * t ) # 產生鋸齒波

plt.plot( t, x )                         # 繪圖
plt.xlabel( 't (second)' )
plt.ylabel( 'Amplitude' )
plt.axis( [ 0, 1, -1.2, 1.2 ] )

plt.show( )

三角波 Triangle Wave

A = 1, f=5Hz, t=0~1s

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

t = np.linspace( 0, 1, 1000 )                 # 定義時間陣列
x = signal.sawtooth( 2 * np.pi * 5 * t, 0.5 ) # 產生三角波

plt.plot( t, x )                              # 繪圖
plt.xlabel( 't (second)' )
plt.ylabel( 'Amplitude' )
plt.axis( [ 0, 1, -1.2, 1.2 ] )

plt.show( )

諧波 Harmonic

兩個以上的弦波組成

\(x(t) = \sum_{k=1}^{N}A_k cos(2 \pi f_k t)\)

\(A_k\) 是第 k 個弦波的振幅，\(f_1\) 稱為基礎頻率 Fundamental Frequency (基頻)， \(f_k\) 是第 k 個弦波的頻率，是 \(f_1\)的整數倍。第一個弦波稱為基礎弦波 Fundamental Sinusoid

ex:

\(x_1(t) = cos(2 \pi 2 t)\)

\(x_2(t) = cos(2 \pi 4 t)\)

諧波 \(x(t) = x_1(t)+x_2(t)\)

import numpy as np
import matplotlib.pyplot as plt

t = np.linspace( 0, 1, 1000, endpoint = False ) # 定義時間陣列

f1 = 2                                          # 定義基礎頻率
x1 = np.cos( 2 * np.pi * f1 * t )               # 產生第1個弦波
x2 = np.cos( 2 * np.pi * 2 * f1 * t )           # 產生第2個弦波
x = x1 + x2                                     # 產生諧波

plt.subplot(131)
plt.plot( t, x1 )
plt.xlabel( 't (second)' )
plt.ylabel( 'Amplitude' )
plt.axis( [ 0, 1, -2, 2 ] )

plt.subplot(132)
plt.plot( t, x2 )
plt.xlabel( 't (second)' )
plt.ylabel( 'Amplitude' )
plt.axis( [ 0, 1, -2, 2 ] )

plt.subplot(133)
plt.plot( t, x )
plt.xlabel( 't (second)' )
plt.ylabel( 'Amplitude' )
plt.axis( [ 0, 1, -2, 2 ] )

plt.show( )

節拍波 beat

兩個不同頻率的弦波，相乘得來

\(x(t) = A cos(2 \pi f_1 t) \cdot cos(2 \pi f_2 t)\) ，A 為振幅， \(f_1\) 是低頻訊號的頻率，\(f_2\)是高頻訊號的頻率

低頻的訊號，構成訊號的 Envelop，高頻訊號，稱為載波訊號 Carrier Signal

ex: \(x(t) = A cos(2 \pi 20 t) \cdot cos(2 \pi 200 t)\)

import numpy as np
import matplotlib.pyplot as plt

t = np.linspace( 0, 0.1, 1000, endpoint = False )   # 定義時間陣列    

f1 = 20                                             # 低頻頻率
f2 = 200                                            # 高頻頻率
x = np.cos( 2 * np.pi * f1 * t ) * np.cos( 2 * np.pi * f2 * t )
envelop1 =  np.cos( 2 * np.pi * f1 * t )            # 包絡
envelop2 = -np.cos( 2 * np.pi * f1 * t )

plt.plot( t, x, '-' )                               # 繪圖
plt.plot( t, envelop1, '--', color = 'b' )
plt.plot( t, envelop2, '--', color = 'b' )
plt.xlabel( 't (second)' )
plt.ylabel( 'Amplitude' )
plt.axis( [ 0, 0.1, -1, 1 ] )

plt.show( )

振幅調變

調幅 Amplitude Modulation (AM) 可定義為

\(y(t)=x(t) \cdot cos(2 \pi f_ct)\) 其中 x(t) 是輸入訊號， \(f_c\) 稱為載波頻率 Carrier Frequency

跟 beat 乘法運算類似。常應用於 AM 收音機、無線電對講機。通常輸入訊號 x(t) 頻率較低，屬於人類的聽力範圍 20Hz ~ 20kHz，載波頻率比較高，例如低頻無線電波，落在 30k ~ 300k Hz

非週期性訊號

淡入與淡出

隨著時間，讓 amplitude 逐漸增加/減少的訊號

ex: \(x(t) = cos(2 \pi 5 t)\) 套用淡入、淡出效果

import numpy as np
import matplotlib.pyplot as plt

t = np.linspace( 0, 1, 1000, endpoint = False ) # 定義時間陣列

x = np.cos( 2 * np.pi * 5 * t )                 # 產生弦波
a = np.linspace( 1, 0, 1000, endpoint = False ) # 產生淡出陣列 
x = x * a                                       # 套用淡出效果 

plt.plot( t, x )                                # 繪圖
plt.xlabel( 't (second)' )
plt.ylabel( 'Amplitude' )
plt.axis( [ 0, 1, -1.2, 1.2 ] )

plt.show( )

剛剛是用線性方式產生淡出陣列，也可以用指數衰減方式產生

\(x(t) = A e^{-αt} cos(2 \pi f t)\)

A 為振幅，f 是頻率，α 為衰減參數，可調整作用時間

啁啾訊號 Chirp

隨著時間，頻率逐漸增加/減少的訊號，分別稱為 Up-Chirp/Down-Chirp

Chirp 經常稱為 Sweep Singal

常應用於 radar 或 sonar，在不同時間點掃描不同的頻率區間/頻帶，藉以偵測目標物

ex: 產生 0~5 秒，由 0Hz 線性增加至 5Hz

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

t = np.linspace( 0, 5, 1000, endpoint = False ) # 定義時間陣列
x = signal.chirp( t, 0, 5, 5, 'linear' )        # 產生啁啾訊號

plt.plot( t, x )                                # 繪圖
plt.xlabel( 't (second)' )
plt.ylabel( 'Amplitude' )
plt.axis( [ 0, 5, -1.5, 1.5 ] )

plt.show( )

wave

需要的套件

pip install numpy
pip install scipy

程式裡面有以下這些波形 wav 的產生方式

• sinusoid() 弦波
• square() 方波
• sawtooth() 鋸齒波
• triangle() 三角波
• harmonic() 諧波
• beat() 節拍波
• fadeout() 淡出
• chirp() 啁啾

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( )

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

        amplitude = 30000,          # 振幅
        frequency = 100,            # 頻率(Hz)
        duration = 3,               # 時間長度(秒)
        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 );

    write_wav_file(wav_params, x);

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

        amplitude = 30000,          # 振幅
        frequency = 100,            # 頻率(Hz)
        duration = 3,               # 時間長度(秒)
        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 * signal.square( 2 * np.pi * wav_params.frequency * t );

    write_wav_file(wav_params, x);

def sawtooth():
    # sawtooth 鋸齒波
    wav_params = WavFileParams(
        filename = "sawtooth.wav",  # 檔案名稱

        amplitude = 30000,          # 振幅
        frequency = 100,            # 頻率(Hz)
        duration = 3,               # 時間長度(秒)
        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 * signal.sawtooth( 2 * np.pi * wav_params.frequency * t )

    write_wav_file(wav_params, x);

def triangle():
    # triangle 三角波
    wav_params = WavFileParams(
        filename = "triangle.wav",  # 檔案名稱

        amplitude = 30000,          # 振幅
        frequency = 100,            # 頻率(Hz)
        duration = 3,               # 時間長度(秒)
        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 * signal.sawtooth( 2 * np.pi * wav_params.frequency * t, 0.5 )

    write_wav_file(wav_params, x);

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

        amplitude = 10000,          # 振幅
        frequency = 100,            # 頻率(Hz)
        duration = 3,               # 時間長度(秒)
        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 );
    x1 = wav_params.amplitude * np.cos( 2 * np.pi * wav_params.frequency * t );
    x2 = wav_params.amplitude * np.cos( 2 * np.pi * (2 * wav_params.frequency) * t );

    x = x1 + x2

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

    write_wav_file(wav_params, x_clip);

def beat():
    # beat 節拍波
    wav_params = WavFileParams(
        filename = "beat.wav",      # 檔案名稱

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

        num_channels = 1,           # 通道數
        sampwidth = 2,              # 樣本寬度
        # num_frames = num_samples, # 音框數 = 樣本數
        comptype = "NONE",          # 壓縮型態
        compname = "not compressed" # 無壓縮
    )
    f1 = 20                     # 低頻頻率(Hz)
    f2 = 200                    # 高頻頻率(Hz)

    t = np.linspace( 0, wav_params.duration, wav_params.num_samples, endpoint = False )
    x = wav_params.amplitude * np.cos( 2 * np.pi * f1 * t ) * np.cos( 2 * np.pi * f2 * t )

    write_wav_file(wav_params, x);

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

        amplitude = 30000,          # 振幅
        frequency = 300,            # 頻率(Hz)
        duration = 3,               # 時間長度(秒)
        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 )
    # 產生淡出陣列, 1 -> 0, num_samples 個值
    a = np.linspace( 1, 0, wav_params.num_samples, endpoint = False )
    x = wav_params.amplitude * a * np.cos( 2 * np.pi * wav_params.frequency * t )

    write_wav_file(wav_params, x);

def chirp():
    # chirp 啁啾, 頻率隨著時間, 逐漸增加或減少
    wav_params = WavFileParams(
        filename = "chirp.wav",     # 檔案名稱

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

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

    f0 = 0                      # 初始頻率(Hz)
    f1 = 1000                   # 終止頻率(Hz)

    t = np.linspace( 0, wav_params.duration, wav_params.num_samples, endpoint = False )
    # 以線性方式遞增頻率
    x = wav_params.amplitude * signal.chirp( t, f0, wav_params.duration, f1, 'linear' )

    write_wav_file(wav_params, x);

def main():
    sinusoid()
    square()
    sawtooth()
    triangle()
    harmonic()
    beat()
    fadeout()
    chirp()

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

References

數位訊號處理：Python程式實作（附範例光碟）（第二版）