第二節 - 訊號處理範例展示

這邊我們試著用 MATLAB 這個程式語言來示範濾波，並且與第二章的理論相互驗證。值得注意的是，一旦將運算搬進計算機，不論時域或頻域，都必須以離散的形式呈現，例如 CTFT 的積分對程式語言來說只是一種極限形式的求和，下一章會繼續探討離散時域下的數學分析 (有些東西現在難以想像)。其實 MATLAB 有很多專門作訊號處理的內建函式，不過為了概念上的清晰，我們這邊先一概不用。

大致步驟如下： (1) 產生範例弦波 (2) 計算 CTFT 分析弦波成分 (3) 觀察頻譜以決定欲留下的頻率成分 (4) 設計濾波函數 (5) 將濾波函數與範例弦波做摺積運算 (6) 獲得過濾後的弦波！ (7) 再次做 CTFT 驗證弦波成分是我們要的 (8) 濾波成果比較

Step 1. 產生範例弦波

我們想觀察的訊號是 x(t)=\begin{cases} A_1\sin(2\pi f_1t)+A_2\sin(2\pi f_2t)&\text{ for }0\le t\le 6 \\ 0&\text{ otherwise,} \end{cases}\quad\text{where}\begin{array}\ (f_1,A_1)=(1,10),\\ (f_2,A_2)=(8,2).\color{white}0\end{array}

翻譯成白話文，這個訊號是由兩種不同頻率的弦波組成，主要頻率 1 Hz 的振幅為 10，而次要頻率 8 Hz 的振幅為 2，但是我們不可能計算一條無限長的訊號，因此只取其中的 6 秒。

我們先用 MATLAB 程式碼產生出 $x(t)$ 。

dt = 0.005; % 時間軸間隔
t_axis = (0:dt:6); % 時間軸，從 0 到 6 秒，每 dt 秒為一間隔
x = 10*sin(2*pi*1*t_axis) + 2*sin(2*pi*8*t_axis); % 波型

% 以下畫圖
figure;
plot(t_axis,x);
xlabel('time (s)');
ylabel('x(t)');
title('Waveform of x(t)');

可觀察到除了有主頻率 1 Hz 的大波外，每一秒內也各有 8 個小波峰與小波谷，代表也隱含 8 Hz 的小波。

Step 2. 計算 CTFT 分析弦波成分

f_axis = (-10:0.005:10); % 頻率軸，從 -10 到 10 Hz，每 0.005 Hz 為一間隔
X = zeros(1,length(f_axis)); % 預先配置好歸零陣列空間，能增進效率
for iFreq = 1:length(f_axis)
  for iTime = 1:length(t_axis)
    % CTFT 公式
    X(iFreq) = X(iFreq) + x(iTime)*exp(-1i*2*pi*f_axis(iFreq)*t_axis(iTime))*dt;
  end
end

mag_X = abs(X); % 計算頻譜量值 (magnitude)

% 以下畫圖
figure;
plot(f_axis, mag_X);
yticks((0:6:36));
ylim([0 36]);
xlabel('Frequency (Hz)');
ylabel('abs(X(f))');
title('Magnitude Spectrum of X(f)');

除了肉眼觀察頻率成分之外，更可靠的方式就是計算 CTFT！結果如上圖，的確是有 1Hz、8Hz 兩種成分。在這邊讀者可以額外思考兩個問題：為什麼 1 Hz、-1 Hz 的量值各為 30，而 8 Hz、-8 Hz 的量值各為 6？又為什麼計算出來的頻譜除了主成分之外，在其他頻率成分也會有少許的值 (sidelobe)？

Step 3. 觀察頻譜以決定欲留下的頻率成分

觀察一下頻譜我們發現，這個複合波的確包含了 1 Hz、8 Hz 兩個頻率成分。由於 8 Hz 的振幅較小，我們姑且假設它是雜訊並且想要加以濾除，因此必須設計一個低通濾波器 (low-pass filter)。所謂的低通就是只讓低頻率波通過的意思，因為我們只需要濾除 8 Hz 的成分，姑且把截止頻率設定在 4 Hz。

Step 4. 設計濾波函數

要設計一個只能讓小於 4 Hz 的波通過的低通濾波器 $h(t)$ ，我們得想辦法讓濾波器頻譜長得像下面這樣：

df = 0.005; % 頻率軸間隔
f_axis = (-5:df:5); % 頻率軸，從 -5 到 5 Hz，每 df 赫茲為一間隔
H = heaviside(f_axis+4) - heaviside(f_axis-4); % 頻譜

% 以下畫圖
figure;
plot(f_axis, H, 'linewidth', 3);
xlabel('Frequency (Hz)');
ylabel('abs(H(f))');
title('Magnitude Spectrum of H(f)');

意思是說在頻域環境下將原本弦波頻譜 $X(f)$ 與濾波器頻譜 $H(f)$ 相乘，便能讓新輸出訊號的頻譜 $Y(f)$ 只保留 1 Hz 此成分。雖然圖片中的頻率軸是從 -10 Hz 到 10 Hz，不過由於大於 4 Hz (小於 -4 Hz) 的值都是零，並不會改變積分值，所以等會兒在計算時頻率軸還是會設為 -5 Hz ~ 5 Hz，以增進效率。我們可以先用 MATLAB 計算一下反轉換 (inverse CTFT，簡稱 ICTFT，把大寫 I 加在前面表示反轉換的寫法後面也會很常出現)，以觀察 $h(t)$ 訊號。這邊要特別注意一下時間單位要跟 $x(t)$ 的一樣，等一下會解釋原因。

t_axis = (-10:0.005:10); % 時間軸，從 -10 到 10 秒，每 0.005 秒為一間隔
h = zeros(1,length(t_axis)); % 預先配置好歸零陣列空間，能增進效率
for iTime = 1:length(t_axis)
  for iFreq = 1:length(f_axis)
    % ICTFT 公式
    h(iTime) = h(iTime) + H(iFreq)*exp(1i*2*pi*f_axis(iFreq)*t_axis(iTime))*df;
  end
end
h = real(h); % 捨棄微小的虛部殘值，這邊是因為我已經知道答案是 sinc 的實訊號才能這樣做
% 另外，從頻域轉回來時域的訊號未必都是實訊號，如果出現複數訊號的話代表此濾波器不符合現實

% 以下畫圖
figure;
plot(t_axis,h);
xlabel('time (s)');
ylabel('h(t)');
title('Waveform of h(t)');

這邊要稍微提醒的是理想的 $h(t)$ 其實是個非時限訊號，電腦上並不能計算，所以說我們要適當地手動近似成時限訊號，一個方便但不太正當的做法是把函數值很小的點直接拿掉，只不過這個函數可以忽略的地方恰好就是時間離 $0$ 很遠的地方，於是這個動作不用特別再做一次了。

Step 5. 將 $x(t)$ 與 $h(t)$ 做摺積運算

回憶一下前面的公式，頻域乘法相當於時域摺積，也就是說我們的輸出訊號 $y(t)=x(t)*h(t)$ 所對應到的頻譜就有上述濾波後的結果。因為電腦計算積分只能用級數和逼近的方式進行，下面示範如何自行實作摺積公式 [2]。

m = length(x); % (0:0.005:6)
n = length(h); % (-10:0.005:10)
x_ext = [x,zeros(1,n-1)]; % 因為摺積會將時間延後，故特地將訊號補零
h_ext = [h,zeros(1,m-1)]; % 因為摺積會將時間延後，故特地將訊號補零
y = zeros(1,n+m-1); % (-10:0.005:16)
for i = 1 : n+m-1
  for j = 1 : m
    if(i-j+1>0)
      y(i) = y(i) + x_ext(j) * h_ext(i-j+1) * dt; % 摺積公式，dt = 0.005
    else
      break % 增進效率
    end
  end
end

Step 6. 獲得過濾後的弦波！

t_axis = (-10:dt:16);
% 以下畫圖
figure;
plot(t_axis,y);
xlim([-2 8]);
xlabel('time (s)');
title('Waveform of y(t)');

摺積後的圖如下：

觀察一下發現 8 Hz 的小波雜訊真的不見了！

Step 7. 再次做 CTFT 驗證弦波成分是我們要的

f_axis = (-10:0.005:10); % 頻率軸，從 -10 到 10 Hz，每 0.005 Hz 為一間隔
Y = zeros(1,length(f_axis)); % 預先配置好歸零陣列空間，能增進效率
for iFreq = 1:length(f_axis)
  for iTime = 1:length(t_axis)
    % CTFT 公式
    Y(iFreq) = Y(iFreq) + y(iTime)*exp(-1i*2*pi*f_axis(iFreq)*t_axis(iTime))*dt;
  end
end

mag_Y = abs(Y); % 量值 (magnitude)

% 以下畫圖
figure;
plot(f_axis, mag_Y);
xlabel('Frequency (Hz)');
ylabel('abs(Y(f))');
title('Magnitude Spectrum of Y(f)');