这里按照题目的基本要求讲解设计流程。

首先通过matlab自带的函数randn函数参数随机数，这里采用基于系统时钟的随机数种子，在matlab中使用的代码语句为：

RandStream.setDefaultStream(RandStream('mt19937ar','seed',sum(100*clock)));

这里使用了RandStream命令实现随机数种子的产生，其中sum(100*clock)为获取系统的时钟，这里如果需要每次产生固定的随机数，只需要将clock修改为一个固定的值就可以了。在完成这个语句之后，执行：

y1 = randn([N,1]);

    为了设计需要，这里还需要对输入的信号进行FFT变换，对于FFT，主要使用MATLAB的自带函数fft进行分析，matlab的代码如下所示：

fy1 = fft(y1,NFFT);

f = fs/2*linspace(0,1,NFFT/2+1);

subplot(322);

plot(f,abs(fy1(1:NFFT/2+1)));

title('Single-Sided Amplitude Spectrum of y(t)')

xlabel('Frequency (Hz)')

ylabel('|Y(f)|')

grid on;

axis([0,NFFT/2,1.2*min(abs(fy1)),1.2*max(abs(fy1))]);

从仿真结果如下所示，随机产生的EMG信号，其频谱图在频域上的各个频率点的都有分布。

function  [EMG_vector60,EMG_vector]=emg_sim(times,power_line_interference);
%times                  ：仿真时间
%power_line_interference：电力线路干扰值
%EMG_vector60           :包含60hz干扰
%EMG_vector             :滤掉60hz干扰


%参数初始化
fs = 2048;    %采样频率
N  = times*fs;%采样点数
%first产生一个随机数向量
RandStream.setDefaultStream(RandStream('mt19937ar','seed',sum(100*clock)));
y1 = randn([N,1]);

figure;
subplot(321);plot(y1,'b');title('an independent, zero-mean, Gaussian random vector');
xlabel('simulation times');
ylabel('amplititude');
grid on;
axis([0,N,1.2*min(y1),1.2*max(y1)]);
NFFT = N; 
fy1 = fft(y1,NFFT);
f = fs/2*linspace(0,1,NFFT/2+1);
subplot(322);
plot(f,abs(fy1(1:NFFT/2+1))); 
title('Single-Sided Amplitude Spectrum of y(t)')
xlabel('Frequency (Hz)')
ylabel('|Y(f)|')
grid on;
axis([0,NFFT/2,1.2*min(abs(fy1)),1.2*max(abs(fy1))]);


%second通过一个巴特沃斯滤波器
%分别为通带截止频率和阻带截止频率 
Wp=21/150;%cut-off:150hz
Ws=40/150;
Rp=1;   %通带最大衰减Rp=3dB;
As=120; %阻带最小衰减As=60              
[n,Wn]=buttord(Wp,Ws, Rp,As);    
[b, a]=butter(n,Wn,'low');                     
%进行滤波
y2=filter(b,a,y1); 
y2=y2*abs(max(y1))/abs(max(y2));
 
subplot(323);plot(y2,'b');title('after filter');
xlabel('simulation times');
ylabel('amplititude');
grid on;
axis([0,N,1.2*min(y2),1.2*max(y2)]);
NFFT = N; 
fy2 = fft(y2,NFFT);
f = fs/2*linspace(0,1,NFFT/2+1);
subplot(324);
plot(f,abs(fy2(1:NFFT/2+1)));hold on;
plot(150,0:1:1000,'r');
title('Single-Sided Amplitude Spectrum of y(t)')
xlabel('Frequency (Hz)')
ylabel('|Y(f)|')
grid on;
axis([0,NFFT/2,1.2*min(abs(fy2)),1.2*max(abs(fy2))]);

EMG_vector=y2;

%加入一个60hz的sin
f_sin=60;
t=1/f_sin/(N/f_sin) : 1/f_sin/(N/f_sin) : N/f_sin/(N/f_sin);
y_sin=power_line_interference*sin(2*pi*f_sin*t);
y3=y2+y_sin';
subplot(325);plot(y3);title('Add sin_600Hz');
grid on;
axis([0,N,1.2*min(y3),1.2*max(y3)]);
NFFT = N; 
fy3 = fft(y3,NFFT);
f = fs/2*linspace(0,1,NFFT/2+1);
subplot(326);plot(f,abs(fy3(1:NFFT/2+1))) 
title('Single-Sided Amplitude Spectrum of y(t)')
xlabel('Frequency (Hz)')
ylabel('|Y(f)|')
grid on;
axis([0,NFFT/2,1.2*min(abs(fy3)),1.2*max(abs(fy3))]);

EMG_vector60=y3;