%問題: 設信號x(t)=sin(2*pi*80*t)+2*sin(2*pi*150*t), 由于某一原因，
%原始信號被白噪聲污染,實際獲得的信號為x2(t)=x+rand(size(x))，
%要求設計一個FIR濾波器恢復出原始信號。
clear all;
close all;
t=0.0:0.001:2.047;
x=sin(2*pi*80.*t);%+2*sin(2*pi*150.*t);
x1=x+randn(size(x));
L=length(x);
n1=1:1024;
n2=1:200;
M=64;

Y=fft(x1);  %噪聲信號的頻譜
E=fft(x);   %原始信號的頻譜
E=abs(E(n1)); 
Y=abs(Y(n1));  %取絕對值

fs=1000;
df=fs/L;       
Wn=[75 80 145 155]/500;
b=fir1(M,Wn);       %hamming window,b是fir濾波器的系數
%freqz(b,1,512);
n=(1:length(b));
stem(n,b);

figure;
z=filter(b,1,x1);   %濾波后的信號
zk=fft(z);          %濾波后信號的頻譜 
zk=abs(zk(n1));
subplot(311);
%plot(n2,x(1536+n2));title('原始信號x');
plot(1:L,x);title('原始信號x');
subplot(312);
plot(1:L,x1);title('在原始信號上加上噪聲信號x');
subplot(313);
plot(1:L,z);title('經過濾波后的信號z');

figure;
subplot(311);
plot(n1*df,E);title('原始信號頻譜');
subplot(312);
plot(n1*df,Y);title('噪聲信號的頻譜');
subplot(313);
plot(n1*df,zk);title('經過濾波后的信號頻譜');