目錄
一、理論基礎
二、核心程式
三、測試結果
一、理論基礎
1.1巴特沃斯濾波器
MATLAB函數“emg_sim（）”：編寫一個MATLAB功能，模拟以2048 Hz采樣的emg信号。函數的輸入應為模拟持續時間（秒）。第二個輸入應用于指定相對于“幹淨”肌電信号RMS值的電力線幹擾幅度。在包含60 Hz幹擾之前，函數的一個輸出應為模拟EMG向量（見下文）。該函數的另一個輸出應該是模拟EMG向量，包括60 Hz幹擾。如果您選擇，可以使用其他參數。要模拟肌電信号，首先使用MATLAB函數“randn（）”建立一個獨立的零均值高斯随機向量。在使用rand（）生成器之前，請確定使用系統時鐘随機地對其進行狀态（也稱為“種子”/“狀态”）。[請參見MATLAB幫助頁中的randn（），以找到一個單行指令。]其次，設計一個截止頻率為150 Hz的二階低通巴特沃斯濾波器。第三，将随機序列通過巴特沃斯濾波器。這樣做可以塑造随機序列，使其看起來像一個持續努力的肌電信号，其頻譜類似于上圖。第四，将适當縮放的60 Hz正弦波添加到信号中（以模拟電力線幹擾）。第五，注意巴特沃斯過濾器将有一個與之相關的啟動瞬态。這種瞬态降低了模拟的保真度。是以，由randn（）建立的初始随機序列的長度必須等于所需的EMG序列長度加上啟動瞬态的長度。在第五步中，“切斷”第四階段序列的啟動瞬态。您的模拟已完成。注意，EMG信号的絕對縮放對于這個問題并不重要。僅EMG信号相對于60Hz幹擾的相對振幅。由于肌電信号是随機的，通常最好将60 Hz幹擾的幅度縮放到模拟“幹淨”肌電的RMS值（低通濾波後）。
使用matlab中自帶的randn函數産生一組随機數，作為EMG信号，然後EMG信号的采樣率為2048hz。這裡随機數産生的随機數種子采用的機遇系統時鐘的随機數種子。系統輸入有兩個，一個是仿真時間，機關為s，一個是幹擾值，輸出有兩個，一個為EMG信号，一格式帶60hz正弦幹擾的EMG信号。系統的設計步驟如下所示：
首先使用randn産生一組随機數，然後設計一個低通的巴特沃斯濾波器，其截止頻率為150hz，将建立的随機數EMG信号輸入到巴特沃斯濾波器。再将濾波得到的信号添加一個60hz的sin信号作為幹擾。進而實作函數一個主要功能。
由于信号隻保留EMG信号中的150hz，而對于150hz之外的頻率信号需要濾除，巴特沃斯(Butterworth)濾波器是一種具有最大平坦幅度響應的低通濾波器,它在通信領域裡已有廣泛應用，在電測中也具有廣泛的用途，可以作檢測信号的濾波器。
巴特沃斯低通濾波器的平方幅度響應為：
buttord和butter來設計巴特沃斯濾波器。buttord函數可在給定濾波器性能的情況下，選擇巴特沃斯濾波器的階數n和截止頻率，進而可利用butter函數設計巴特沃斯濾波器的傳遞函數。
1.2 陷波器
陷波器是一種特殊的帶阻濾波器,其阻帶在理想情況下隻有一個頻率點，是以也被稱為點阻濾波器。這種濾波器主要用于消除某個特定頻率的幹擾，由于陷波器頻率特性的特殊性，它除了可采用雙線性變換進行設計外，還可以采用所謂零極點配置的方法進行設計。
一個理想的陷波濾波器的頻率特性要在消除的信号頻率點處，其值等于零；而在其他頻率點處，其值等于1。由于數字濾波器的頻率特性就是其機關沖激響應在機關圓上的Z變換，是以隻需要在機關圓上相應于所需帶阻濾波器阻帶位置的頻率處設定零點，就可以使濾波器的頻率特性在所需阻帶頻率處為零。但是僅僅進行零點設定隻考慮到了濾波器的阻帶特性。為了得到非常陡峭的過渡帶和常數幅度的通帶特性，必須在Z平面上為每一個零點再配置一個相應的極點。Z平面機關圓附近的零點會在濾波器幅頻特性的相應頻率處産生陷落，零點離機關圓越近，陷落越深；而Z平面機關圓附近的極點會在濾波器幅頻特性的相應頻率處産生凸峰，極點離機關圓越近，凸峰越高。是以在完成了零點的配置後，為了抵消零點引起的陷落對濾波器通帶範圍内幅頻特性的影響，還需要再配置相應的極點，由于濾波器穩定性的要求，極點必需配置在機關圓内，顯然極點離機關圓越近則極點對零點的抵消作用越明顯，得到的濾波器的阻帶就越窄，過渡帶就越陡峭。
二、核心程式
clc;
clear;
close all;
times = 1;
power_line_interference = 1;
order = 1;
%function 1：generator EMG
[EMG_vector60,EMG_vector] = emg_sim(times,power_line_interference);
%function 2：notch_filter
EMG_notch_Filter = emg_notch60(EMG_vector60,order);
%function 3
EMG_notch = emg_freq_null(EMG_vector60);
%function 4
RMS = emg_compare(EMG_vector,EMG_notch); function EMG_notch_Filter = emg_notch60(EMG_vector60,order);
%EMG_60hz:輸入60hz的MEG信号
%order :notch濾波器的階數
fs = 2048;
N = length(EMG_vector60);
Wp = [1 150]/1000;
Ws = [59 61]/1000;
Rp = 1;
Rs = 15;
[n,Wn] = buttord(Wp,Ws,Rp,Rs);
[b2,a2] = butter(order,Wn,'stop');
figure
freqz(b2,a2,1000,2000)
title('n=2 Butterworth Bandstop Filter')
EMG_notch_Filter=filter(b2,a2,EMG_vector60);
figure
%加入一個60hz的sin
subplot(211);plot(EMG_vector60);title('EMG+sin');grid on;
axis([0,N,1.2*min(EMG_vector60),1.2*max(EMG_vector60)]);
subplot(212);plot(EMG_notch_Filter);title('EMG-sin');grid on;
axis([0,N,1.2*min(EMG_notch_Filter),1.2*max(EMG_notch_Filter)]); 三、測試結果
