clear all;                    	     				   	%清寄存器值
clf;                              				   	%清屏
N=256;                            			 	   	%數據點數
fs=10000;                            			        	%采樣頻率
dt=1/fs;                                		    		%采樣時間間隔
tt=0.:dt:dt*(N-1);
dis=(rand(size(tt)));                              		    	%隨機噪聲
for k=1:N;                    
   f1=3000;                         			 		%信號1頻率
   f2=1000;                         			 		%信號2頻率
   y(k)=sin(2*pi*f1*k*dt)+sin(2*pi*f2*k*dt)+dis;  %產生合成信號
end
lp=3500;                                   	 		%截止頻率
wn1=2*lp/fs;                              	 			%函數的參數
[z1,p1,k1] = CHEBY1(3,0.5,wn1);            	 			%濾波器的極零點表示
[b1,a1] = CHEBY1(3,0.5,wn1);              	 			%濾波器的傳遞函數表示
b1=b1/(8*1.0711);                    		 			%將參數按比例縮小
a1=a1/(8*1.0711);                             		%為DSP程序做準備，MATLAB中不需要

yy1=filter(b1,a1,y);                   	      			%濾波

y=fft(y,N);                         	      			%將原始信號做FFT變換
pyy=y.*conj(y);                     	      			%做功率譜分析
f=(1:(N/2－1));
for i=1:N/2-1;
   f(i)=f(i)*fs/N;
end
figure(1);

plot(f,pyy(2:N/2))                              			%畫出原始信號的頻譜圖

y=fft(yy1,N);                      	       			%將濾波后數據做功率譜分析
pyy=y.*conj(y);
f=(1:(N/2－1));
for i=1:N/2－1;
   f(i)=f(i)*fs/N;
end
figure(2);
plot(f,pyy(2:N/2));                               			%畫出濾波后信號的頻譜圖