Matlab應用：數字濾波器 貫穿時域、頻域以一階電路 RC Low pass為例 (二)

本篇進行訊號延遲的說明以及轉換為 C 語言的教學。

概述訊號(相位)延遲是什麼

承上篇，訊號改用單一頻率 50Hz 的輸入訊號，一樣經由一階濾波器截止頻率 50 Hz

圖(1) 為，單一頻率 50Hz 的輸入訊號(振幅為10)  【MATLAB】

圖(2) 為，單一頻率 50Hz 的輸入訊號(振幅為10)，但經過了 cut-off frequency = 50 Hz 的一階低通濾波 【MATLAB】

圖(3) 為 圖(1) 及 圖 (2) 比較圖

圖(4) 為 cut-off frequency = 50 Hz 的一階低通濾波 之波德圖

由自動控制理論中的 Bode plot 我們可以得知，當輸入訊號為 50Hz，經過一個一階低通濾波器且截止頻率亦為 50Hz 時，兩個頻率正好會對應到 -3dB 點，對於振幅而言會削弱至原本的 0.707 倍；對於相位而言會產生 45 度的相位延遲(見圖4)。

振幅的部分很好理解，由圖(3)可以看出 50Hz 頻率對應到的振幅由 10 降到了 0.707 ，表示輸入訊號 50Hz 的部分受到了抑制(以此類推，上篇中的 700Hz 被削弱到幾乎消失，這也是"低通"濾波器名稱的由來，"低的頻率導通" "高的頻率被濾除")，但相位這個部分就需要稍微理解一下。

在 MCU 實現的過程，一般而言我們都會將程式利用中斷檔反覆執行，而此時濾波器造成的訊號延遲就有可能會導致輸出不如我們預期，如圖(3)中的上圖，我們可以得知輸入訊號與輸出訊號峰值時間差為0.0024秒(約為0.0025秒)，而我們又知道整個訊號是由 50Hz 的單一頻率所組成，50Hz 的訊號其週期為 1/50 = 0.02秒，由這兩者的關係我們可以知道時間差大略為整個波型的 0.0025/0.02 = 0.125 ，這個數字代表時間延遲佔了整個週期的 12.5% 。我們將其轉換為相角的角度來看，整個波型一週期為 360 度，乘上 0.125 為 45 度，所謂的相位延遲 45 度就是這樣得來的。

這樣的延遲對於 MCU 實現的過程中會有怎樣的影響? 由於我們在撰寫程式時，其演算法會希望是實時輸出(也就是當你程式一執行的當下，你希望馬上輸出響應)，但實際上這是不可能發生的，光是 ADC 就需要花費數個 us 甚至是 ms 等級的時間，但你的硬體還再持續的作動(例如馬達還在持續的旋轉)，此時參數都已經變動，但因為你還在進行此輪的演算，故有可能會造成響應不如預期甚至系統崩潰的可能性(甚至還需要針對延遲時間進行額外的補償)。所以延遲能避免就進量避免，演算法能簡單就盡量簡單。

濾波器以 C 語言實現

延遲的部分很粗淺的講大致上就是這樣，在 Matlab之中我們已經知道該如何使用 filter，但實際使用中我們往往是需要使用 MCU 撰寫 C 語言來實現 filter 的功能，接下來我們要學如何將 Z-domain 轉成差分方程後用 C 語言實現。

首先我們從 Matlab help 中的 filter 函數擷取了以下的圖片，是關於 Moving average 的部分( Moving average 其實大家不陌生，在股市市場就很常見，也就是大家常看到的五日線、十日線等等)

上圖中由上到下，大致上是一開始用 linspace 建立了一個給 x index使用的陣列，而 x 則是原訊號。後續給出了 Moving average 的數學式，以中文語意來說 " 當下的輸出為前幾個輸入的組合並且平均 "

其中的 windowSize = 5 事實上就是代表股市的五日線啦！其推導過程如下圖(由 Z-domain直接轉至時域)，在 Z-domain中，Z的負次方在時域代表的就是延時。

由最後的 y(n) = .... 我們嘗試利用 help 內的資料，不使用 filter 這個 function 來做 moving average ，而是直接編寫。

clc
clear all
close all

t = linspace(-pi,pi,100);
x = sin(t) + 0.25*rand(size(t));
windowSize = 5;
b = (1/windowSize)*ones(1,windowSize)
a = 1;
y = filter(b,a,x);
plot(t,x)
hold on
plot(t,y)
grid on
legend('Input Data','Filtered Data','Location','NorthWest')
title('Plot of Input and Filtered Data')

%%
M = 5
YY(1) = x(1)/M
YY(2) = (x(1)+x(2)) /M
YY(3) = (x(1)+x(2)+x(3))/M
YY(4) = (x(1)+x(2)+x(3)+x(4))/M
for i=5:length(x)
    Y(i) = x(i) + x(i-1) + x(i-2) + x(i-3) + x(i-4)
    YY(i) = Y(i)/M
end
figure
plot(t,x,t,y,t,YY,'-r')

上圖中 data1 是原訊號；data2是使用 filter function ，對原訊號濾波後得到的訊號；data3是自己編寫 moving average 對原訊號濾波後的訊號。 可看出 data2、data3 是重疊的，表示這樣的寫法是可行的。

使用 C 語言撰寫一階低通濾波

首先，建立一個 input signal 振幅為 10 ，頻率為 50 Hz

接著建立一階低通濾波器，並繪製出其頻率響應( Bode plot)

直接使用 MATLAB 裡面的 filter function 進行濾波，得到以下

接著我們對轉移函式使用 MATLAB 裡面的 c2d 得到 z-domain 的差分方程，得到該差分方程後就可以用來寫成 C 語言的形式了

比較使用 MATLAB filter function 濾波過後的訊號，以及使用 c2d 轉換得到的 Z 差分方程，對原始訊號進行濾波後的訊號。可以發現濾波後的訊號非常相似(但兩者之間存在一些相位差，源自於離散以及連續域的差異，這部分不在本文中敘述)

最後，比較兩個方法的濾波器所得到的輸出訊號

可以看出，由 filter function 濾波過後的訊號(左)，以及自己 coding 出來的訊號(右)看不出什麼差異，頻譜分析也一致，驗證了這段代碼的可行性，由於沒有使用 Matlab 內部的函數(filter)來進行濾波，所以這部分的代碼就可以很容易的經由一點點小修改，套用到 C 語言或是 Arduino 程式之中了！

整篇文的關鍵在於 要如何將轉移函式變成 C 語言以利我們在 MCU 中實現 ?
答案就是 

1. 先弄出濾波器的轉移函式
2. 使用 MATLAB 中的 c2d 函式 (第55行)
3. 將第2點得到的 Z 差分方程(sysd1)，移項後就可以得到輸出入之間的關係
4. 將第3點得到的輸出入關係寫成 C 語言(112行)

這兩篇從一開始介紹低通濾波器，到最後的 C 語言實現，無疑是為了後續的 Arduino 教學鋪路，希望還不知道如何實現一階 Low pass 數位濾波的同好可以從這得到一些有用的東西，那就不枉費我弄了快一整天阿 ~

而有些人可能會問，濾波器的型式有很多種，怎麼只有教低通濾波器呢?
實際上 Matlab 有現成的函式可以呼叫，例如 Butterworth 的 butter、Chebyshev IIR可以利用 cheb1ord 定階，再利用 cheby1 就可以直接求出 Z-domain 的分母分子項，後續的大家就可以舉一反三了。

除了有現成的函式可以呼叫外，還有 Toolbox可以直接圖形化的進行設計，待之後有空閒時間再補上。


最後，附上程式碼

% http://wyj-learning.blogspot.tw/2017/08/matlab-rc-low-pass_30.html
clc
clear all
close all

% 建立 input signal ( 50Hz sin函數)
fs = 2048;
Ts = 1/fs;
N = 1024;
n = [0:1:N-1];
freq_resolution = fs/N; % 解析度太低-> N上升
t=n/fs;

f1 = 50;
f2 = 0;
input =  10*sin(2*pi*f1*t) + 1*sin(2*pi*f2*t);
figure
subplot(2,1,1);
plot(t,input)
title('原始訊號(時域)')
xlabel('t');
axis tight

% 原始訊號頻譜分析
NFFT = N;
Y = fft(input,NFFT); % Spectrum
mag=abs(Y);
subplot(2,1,2);
f=n*fs/NFFT;   % n*fs/NFFT，若要提高fft頻率範圍僅有提升fs(in this case)
plot(f(1:NFFT/2),mag(1:NFFT/2)*2/NFFT,'-b'); 
grid on
xlabel('Frequency(Hz)');
ylabel('Magnitude');
title('FFT of original-signal');
axis([-inf,inf,0,15])

%% 1st-order low pass ( cut-off freq. = 50Hz  )
RC = 1/(2*pi*50);

num = [1]; % 分子多項式係數
den = [RC 1]; % 分母多項式係數
Hs_1 = tf (num,den)

% 極零點圖
%pzmap(Hs_1)

% Bode plot
figure
bode(Hs_1)
title('Bode plot of Filter (cut-off frequency = 50Hz)')
grid on

%% 濾波器離散化
sample_time = Ts;
sysd1 = c2d(Hs_1, sample_time) % 轉為 Z-domain 方法默認為zoh法、濾波器大多採用'imp'、控制器大多採用'tustin'
%opt = c2dOptions('Method', 'imp');
%sysd1 = c2d(Hs_1, sample_time, opt)
sysd1_num = sysd1.num{1}; % 分子
sysd1_den = sysd1.den{1}; % 分母

% 查看濾波器頻率響應
figure
[A,W] = freqz(sysd1_num, sysd1_den);
plot(W/pi*fs/2, 10*log(abs(A)))
title('濾波器頻率響應(Z-domain)')
ylabel('Amplitude(dB)');
xlabel('frequency(Hz)');

% 濾波(使用filter function)
filter_input = filter(sysd1_num,sysd1_den,input);

% input t域
figure
subplot(2,1,1);
plot(t,filter_input);
xlabel('t(sec)');
title('濾波後的訊號(時域)')
axis tight

% 頻譜分析
NFFT = N
filter_Y = fft(filter_input,NFFT); % Spectrum, NFFT = 1024
filter_mag=abs(filter_Y);
subplot(2,1,2);
f=(0:NFFT-1)*fs/NFFT;   
plot(f(1:NFFT/2),filter_mag(1:NFFT/2)*2/NFFT,'-b'); 
grid on
xlabel('Frequency(Hz)');
ylabel('Magnitude');
title('FFT of filtered-signal');
axis([-inf,inf,0,15])

%% 結果比較
figure
subplot(2,1,1);
plot(t,input,'-r',t,filter_input,'-b')
title('濾波前(紅色)&濾波後(藍色)的訊號比較(時域)')
xlabel('t');


subplot(2,1,2);
plot(f(1:NFFT/2),mag(1:NFFT/2)*2/NFFT,'-r',f(1:NFFT/2),filter_mag(1:NFFT/2)*2/NFFT,'-b'); 
grid on
xlabel('Frequency(Hz)');
ylabel('Magnitude');
title('FFT of original & filtered-signal');
axis([-inf,inf,0,15])

%% 自行撰寫一階低通濾波器(不用filter function)
my_filter_input(1) = input(1);
for i=2:length(input)
    my_filter_input(i) = 0.8578*my_filter_input(i-1) + 0.1422*input(i);
end
figure
plot(t,input,t,filter_input,t,my_filter_input,'r')
legend('input','filtered-signal(with filter function)','filtered-signal(without filter function)')

figure
subplot(2,2,1);
plot(t,filter_input)
title('filtered-signal(with filter function)');
subplot(2,2,2);
plot(t,my_filter_input,'-r')
title('filtered-signal(without filter function)');

subplot(2,2,3);
NFFT = N
filter_Y = fft(filter_input,NFFT); % Spectrum, NFFT = 1024
filter_mag=abs(filter_Y);
f=(0:NFFT-1)*fs/NFFT;   
plot(f(1:NFFT/2),filter_mag(1:NFFT/2)*2/NFFT,'-b'); 
grid on
xlabel('Frequency(Hz)');
ylabel('Magnitude');
title('FFT of filtered-signal(with filter function)');
axis([-inf,inf,0,15])

subplot(2,2,4);
NFFT = N
my_filter_Y = fft(my_filter_input,NFFT); % Spectrum, NFFT = 1024
my_filter_mag=abs(my_filter_Y);
f=(0:NFFT-1)*fs/NFFT;   
plot(f(1:NFFT/2),my_filter_mag(1:NFFT/2)*2/NFFT,'-r'); 
grid on
xlabel('Frequency(Hz)');
ylabel('Magnitude');
title('FFT of filtered-signal(without filter function)');
axis([-inf,inf,0,15])

補充：
關於 moving average 的頻率響應可見此連結
moving average 更多相關實現可見此連結