[signal, freq]=wavread('orig1.wav'); %read orig1.wav file and return 2 parameters
n=length(signal);

figure(1);
plot((1:n)/freq,signal);
title('original signal in the time domain');
xlabel('second');
grid;

figure(2);
signal10=signal-mean(signal);
fsignal=fft(signal10);
plot((1:n/2)/n*freq,abs(fsignal(1:n/2)));

title('original signal in frequency domain');
xlabel('Hz');
ylabel('magnitude');
grid;

figure(3);
plot((1:n/2)/n*freq,unwrap(angle(fsignal(1:n/2))));
title('original signal in frequency domain');
xlabel('Hz');
ylabel('phase');
grid;

sound(signal,freq); %playing sound

Plotting the noisy signal for time , frequency (magnitude and phase) domain.

[signal, freq]=wavread('noisy1.wav'); %read noisy.wav file and return 2 parameters
n=length(signal);

figure(4);
plot((1:n)/freq,signal);
title('noisy signal in the time domain');
xlabel('second');
grid;

figure(5);
signal10=signal-mean(signal);
fsignal=fft(signal10);
plot((1:n/2)/n*freq,abs(fsignal(1:n/2)));

title('noisy signal in frequency domain');
xlabel('Hz');
ylabel('magnitude');
grid;

figure(6);
plot((1:n/2)/n*freq,unwrap(angle(fsignal(1:n/2))));
title('noisy signal in frequency domain');
xlabel('Hz');
ylabel('phase');
grid;

sound(signal,freq); %playing sound

Filterring process for experiment Part 1 using the following function ( fir1, freqz , filter ) , the calculation of MSE and SNR values and the rest of procedure in part 1

[signal2, freq]=wavread('orig1.wav');
[signal1, freq]=wavread('noisy1.wav');
a=signal1;

fs=3417;
fp=243;
Wn=(fs+fp)/freq; %calculate the normalized cutoff frequency
N=48;

b=fir1(N,Wn);
figure(7);
stem(0:length(b)-1,b);%stem the impulse response
title('Impluse response');
xlabel('samples');

figure(8);
freqz(b,1,N);
y=filter(b,1,a);
M=length(y);

sound(y,freq); %play the filtered signal sound
wavwrite(y,'partone'); %save the filtered signal as partone.wav file

MSE = sum((y(1:M)-signal2(1:M)).^2)/M %calculation of MSE(mean squared error)
NSR = 10*log((sum(signal2(1:M).^2)/M)/MSE) %calculation of SNR(signal to noise ratio)
%calculation for SNR for the noisy signal before filtering
MSE = sum((signal1(1:M)-signal2(1:M)).^2)/M ;
NSR = 10*log((sum(signal2(1:M).^2)/M)/MSE)

Ft= 8000;
%DSP chips coefficient, n
n = 4;
%Quantization program equation
Bk = (floor(b.*(2.^(n-1))+0.5))/(2.^(n-1));
%plot the frequency using quantized coefficients
figure;
freqz(Bk,1,N,Ft)
title('4 bit design filter');

%DSP chips coefficient, n
n = 8;
%Quantization program equation
Bk = (floor(b.*(2.^(n-1))+0.5))/(2.^(n-1));
%plot the frequency using quantized coefficients
figure;
freqz(Bk,1,N,Ft)
title('8 bit design filter');

%DSP chips coefficient, n
n = 12;
%Quantization program equation
Bk = (floor(b.*(2.^(n-1))+0.5))/(2.^(n-1));
%plot the frequency using quantized coefficients
figure;
freqz(Bk,1,N,Ft)
title('12 bit design filter');

%DSP chips coefficient, n
n = 16;
%Quantization program
Bkz = (floor(b.*(2.^(n-1))+0.5))/(2.^(n-1));
%plot the frequency using quantized coefficients
figure;
freqz(Bkz,1,N,Ft)
title('16 bit design filter');

Plotting the filtered signal for time , frequency (magnitude and phase) domain.

[signal, freq]=wavread('partone.wav'); %read the partone.wav file which is the filtered signal
n=length(signal);

figure(9);
plot((1:n)/freq,signal);
title('filtered signal in the time domain');
xlabel('second');
grid;

figure(10);
signal10=signal-mean(signal);
fsignal=fft(signal10);
plot((1:n/2)/n*freq,abs(fsignal(1:n/2)));

title('filtered signal in frequency domain');
xlabel('Hz');
ylabel('magnitude');
grid;

figure(11);
plot((1:n/2)/n*freq,unwrap(angle(fsignal(1:n/2))));
title('filtered signal in frequency domain');
xlabel('Hz');
ylabel('phase');
grid;