%----------------- Filter Design -----------------
%
% Technical University of Sofia - MUSSPI lab course
% Created by Viktor Stoev
% November 2016
%----------------------------------------------------


close all;
clear all;
tic
%% Butterworth, Cheby1, Cheby2 and elliptic prototypes

n = 100;  % filter order 
Rp = 1; % ripple magnitude in dB
[z, p, k] = buttap(n); % try [z, p, k] = cheb1ap(n, Rp); %,  cheb2ap(5, 20);, ellipap(5, 0.5, 20);
figure(1);
plot(p, 'x'); %plot poles
hold on;
plot(z, 'o','linewidth',2); %plot zeros
title('Pole-zero plot');
axis equal;
axis([-1.5 1.5 -1.5 1.5]);

w = 0:0.01:5; % frequency range for frequency response
f1 = w/(2*pi);
[b, a] = zp2tf(z, p, k);
h = freqs(b, a, w);
figure(2) 
plot(f1, abs(h),'LineWidth',2);
grid on
ylabel('Magnitude');
xlabel('Frequency,[Hz]');

figure(3)
plot(f1, unwrap(angle(h))); %use unwrap function
grid on
ylabel('Phase(rad)');
xlabel('Frequency(rad/s)');
title('Phase change');


toc
pause
clear all;

%% Generate real filters based on butter, elliptic, etc. protoypes - use butter, cheby1, cheby2, ellip ;
%  Apply the above filter to this sample speech signal

% hh = input('Close figures [y/n]: ','s');
% 
% while(hh ~= 'y' && hh ~= 'n')
%     hh = input('Close figures [y/n]: ','s');
% end
% 
% if(hh == 'y')
%     close all;
% end

close all;

[y,Fs,nbits]=wavread('speech_dft.wav');
figure;                                  
plot(y);
xlabel('Time');
ylabel('Signal magnitude');
title('Speech signal');
 soundsc(y,Fs,nbits);
 pause

f = Fs*((0:(size(y,1))/2)/(size(y,1)));

% design filter to add noise
a1=[0.84,0.26,-0.31,0.13,0.18,-0.29,0.18];
b1=1;
%apply filter
zz = filter(b1,a1,y);
soundsc(zz,Fs,nbits);
pause

FT = abs(fft(y).^2);
FT2 = abs(fft(zz).^2);

figure;
subplot(2,1,1)
plot(f,FT(1:end/2+1,1));
title('Original signal')
ylabel('Magnitude');
xlabel('Frequency.[Hz]');
subplot(2,1,2)
plot(f,FT2(1:end/2+1,1));
ylabel('Magnitude');
xlabel('Frequency,[Hz]');


% filter design to remove noise
n = 6;        % filter order  
Wn = 4000/Fs; % cut-off freq
[b,a] = butter(n,Wn);
fvtool(b,a);

zz_filtered = filter(b,a,zz);
FT3 = abs(fft(zz_filtered).^2);

figure;
subplot(3,1,1)
plot(f,FT(1:end/2+1,1));
title('Original signal')
ylabel('Magnitude');
xlabel('Frequency.[Hz]');
subplot(3,1,2)
plot(f,FT2(1:end/2+1,1));
ylabel('Magnitude');
xlabel('Frequency,[Hz]');
subplot(3,1,3)
plot(f,FT3(1:end/2+1,1));
ylabel('Magnitude');
xlabel('Frequency,[Hz]');
break
pause
soundsc(zz_filtered,Fs,nbits);
