%程序中設定采樣點數為256個，采樣頻率為10000hz，
%輸入信號為隨機噪聲和兩個正弦的合成信號
%設定的濾波器的截止頻率為3500hz
%通過挈比雪膚濾波器，運行程序，比較濾波前后的頻域波形
%濾波去除了高頻信號
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做準備

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));%畫出濾波后信號的頻譜圖