[a,fs]=wavread('現在語音處理技術.wav');%讀入寧波大學語音
a=a(:,1);
FS=length(a);
f=0:fs/FS:(FS-1)*fs/FS;

%進行FFT變換
x=fft(a);
figure(1)
plot(f,abs(x));title('原始信號a的fft圖');
axis([0 45000 0 1000]);
xlabel('f');
ylabel('abs(x)');
grid on;




%設置FIR濾波器
format long
fedge=input('Band edges in Hz=');
mval=input('Desired magnitude values in each band=');
dev=input('Desired ripple in each band=');
FT=input('Sampling frequency in Hz=');
[N,fpts,mag,wt]=remezord(fedge,mval,dev,FT);
b=remez(N,fpts,mag,wt);
disp('FIR Filter Coefficients');
disp(b)
[h,w]=freqz(b,1,256);


figure(2)
subplot(2,1,1);
plot(w/pi,20*log10(abs(h)));title('幅頻特性'); %畫出幅值
xlabel('\omega/\pi');
ylabel('Gain,dB');
grid on;

subplot(2,1,2)
pha=angle(h); %得到相位
plot(w/pi,pha);title('相頻特性');
xlabel('\omega/\pi');
ylabel('phase(Rad)');



%原始語音信號圖
ts=1/fs;
ta=(length(a)-1)/fs;
t=0:ts:ta;
figure(3)
subplot(3,1,1);
plot(t,a);title('原始信號a');
xlabel('t');
ylabel('a');
grid on;

%添加隨機噪聲
noise_mu = 0; 
noise_var =0.005; 
a0= randn(size(a)) .* sqrt(noise_var) + noise_mu;
a1=a+a0;

%在原始信號上加隨機噪聲的信號圖
subplot(3,1,2);
plot(t,a1);title('在原始信號上加上噪聲信號a1');
xlabel('t');
ylabel('a1');
grid on;


%求卷積
a11=conv(b,a1);
 
%經過濾波器后的信號圖
M=(length(a11)-1)/fs;
m=0:ts:M;
subplot(3,1,3);
plot(m,a11);title('經過濾波器后的信號a11');
xlabel('m');
ylabel('a11');
grid on;