% Digital Audio Course 
% © 2014 Prof. Mario Malcangi

% Finite Impulse Response (FIR) Filter
%
% y(n) = h(0)x(n)+h(1)x(n-1)+h(2)x(n-1)+...
% h(n) = Filter impulse response
% H(k) = Filter frequency response
% 

clear;
close all;

Fs=8000;    % Sampling rate
Ts=1/Fs;    % Sampling period

% --- Compute the impulse response h(n) from frequency response H(k) ------

W = 0.1;    % Impulse response duration (s)
K= 0.1/Ts;  % Samples in the window 
F0 = 1/W;   % Fundamental frequency (Hz)

% Design the frequency response

Fc1 = 0;             % Cutoff frequency
Kc1 = (Fc1/F0)+1;    % Cutoff frequency index

Fc2 = Fs/2;      % Cutoff frequency
Kc2 = Fc2/F0;    % Cutoff frequency index

for k=1:Kc1
    H(k)=0;     % Stop band
end

for k=Kc1+1:Kc2
    H(k)=1;     % Pass band
end

for k=Kc2+1:K/2
    H(k)=0;     % Stop band
end

% ------ h(n)= IDFT(H(k)) -------------------------------------------------

% Set positive and negative frequencies

for k=2:K/2
    Xr(k)= 0.5*H(k);        % Positive frequencies (half amplitude)
    Xi(k)= 0;
    Xr(K-k+2)= 0.5*H(k);   % Negative frequencies (half amplitude)
    Xi(K-k+2)=0;
end

% Start the IDFT (using sin and cos)---------------------------------------
N=K;        
for n=0:N-1      % n is the time index t=n*Ts   
    x(n+1)=0;
    for (k=0:K-1)    % k is the frequency index (f=k*F0)
        s(n+1)=sin(2*pi*F0*k*n*Ts);  % n is the time index (t=n*Ts)
        c(n+1)=cos(2*pi*F0*k*n*Ts);  % n is the time index (t=n*Ts)         
        x(n+1)=x(n+1)+Xi(k+1)*s(n+1)+Xr(k+1)*c(n+1);      % sin table lookup
    end
end

Fscale=0:F0:(Fs/2)-F0;
subplot(4,1,1);
plot(Fscale,H);
legend('Frequency response');

% -------- Build h(n) -----------------------------------------------------
m=1;
for n=N/2+1:N   % Put the second half response in the first half
    h(m)=x(n);
    m=m+1;
end
for n=1:N/2
    h(m)=x(n);  % Put the first half response in the second half
    m=m+1;
end
h=h*(2/N);

Tscale=0:Ts:(length(h)*Ts)-Ts;
subplot(4,1,2);
plot(Tscale,h);
legend('Impulse response');

% -------------------------------------------------------------------------
% Pattern to test the filter frequency response
% Sequence of pure tones

D=100;      % Pattern duration (samples)
f0=100;     % Step frequency

f=0;
for k=1:K
    x(k)=0;
end
k=K+1;
while f<Fs/2 % Stop when frequency is higher than half of sampling rate   
    for n=1:D  % Generate a D-length segment of the pattern
        x(k)=sin(2*pi*n*f/Fs);
        k=k+1;
    end
    f=f+f0; % Step the frequency
end


% ----------- FIR filter ---- ---------------------------------------------

N=length(x);
M=K;
n=K+1;  
while n <= N   % Run the filter for all input samples
    y(n)=0;
    for m=1:M
        y(n)=y(n)+h(m)*x(n-m);   % FIR kernel
    end
    n=n+1;
end

% -------------------------------------------------------------------------

% Show how input is filtered
% Check that cutoff frequency is down by -3dB (0.707)
% Check that each octave after cutoff frequency is down by -6dB (0.5)
% Check that all frequencies befor cutoff frequency are always not down 
% more then -3dB

Tscale=0:Ts:(length(x)+100*Ts)-Ts;
subplot(4,1,3)
plot(x(M+1+D:N-1));
legend('Filter Input');
Tscale=0:Ts:(length(y(M+1+500:N-1))*Ts)-Ts;
subplot(4,1,4)
plot(y(M+1+D:N-1));
legend('Filter Output');


  
