Plotting the Graphs
[signal,freq]= wavread('orig1.wav');
n = length (signal);

[noisy,freq]= wavread('noisy1.wav');
m = length (noisy);

Fp=1200;
Fs=1800;
Ft=8000;
Fc=((Fp+Fs)/2);
Ws=2*Fs/Ft;
Wp=2*Fp/Ft;
Wn=(2*Fc/Ft);

[N,Wn]=buttord(Wp,Ws,0.03,40);
disp(N);
[B,A]=butter(N,Wn);
freqz(B,A,N,8000);
y=filter(B,A,noisy);
wavwrite(y,'butfil.wav');

figure;
stem(B);
title('filter impulse response');
xlabel('index,n');
ylabel('impulse response');
grid;
 

Filter frequency response
figure;
[H,W]=freqz(B,A);
magH=abs(H);
db=20*log10(magH);
plot(W,db);
title('Frequency Response of the lowpass filter');
xlabel('normalised frequency,Hz');
ylabel('Gain,dB');
grid;

figure;
plot(W,unwrap(angle(H)));
title('Phase Frequency Response');
xlabel('normalised frequency,Hz');
ylabel('phase');
grid;
 

%MSE of the original signal and the filtered signal
MSE= sum((y-signal).^2)/n;
disp(MSE);

%SNR of the original signal and the filtered signal
SNR=10.*log10(1/n.*(sum(signal.^2))./MSE);
disp(SNR);

%MSE of the original signal and the noisy signal
MSE=sum((noisy-signal).^2)/n;
disp(MSE);

%SNR of the original signal and the noisy signal
SNR=10.*log10(1/n.*(sum(signal.^2))./MSE);
disp(SNR);
 

%comparison between the filtered and original signal in time domain
figure;
subplot(211);
plot ((1:n)/freq,signal);
title('Original signal in time domain');
xlabel ('second');
grid;

subplot(212);
plot ((1:n)/freq,y);
title('Filtered signal in time domain');
xlabel ('second');
grid;
 

%comparison between the filtered and original signal in frequency domain (magnitude)
figure;
signal0=signal -mean (signal);
fsignal=fft(signal0);
subplot(211);
plot ((1:n/2)/n*freq, abs(fsignal(1:n/2)));
title('Original signal in frequency domain');
xlabel ('Hz');
ylabel('Magnitude');
grid;

signal1=y -mean (y);
fsignal1=fft(signal1);
subplot(212);
plot ((1:n/2)/n*freq, abs(fsignal1(1:n/2)));
title('Filtered signal in frequency domain');
xlabel ('Hz');
ylabel('Magnitude');
grid;

%comparison between the filtered and original signal in frequency domain (phase)
figure;
subplot(211);
plot((1:n/2)/n*freq, unwrap(angle(fsignal(1:n/2))));
title ('Original signal in frequency domain');
xlabel ('Hz');
ylabel ('Phase');
grid;

subplot(212);
plot((1:n/2)/n*freq, unwrap(angle(fsignal1(1:n/2))));
title ('Filtered signal in frequency domain');
xlabel ('Hz');
ylabel ('Phase');
grid;