% Programmazione e architetture DSP
% Prof. Mario Malcangi


% Modellazione e simulazione dei modelli computazionali dei sistemi lineari
% (Filtraggio IIR)

% Si implementa un IIR (RC) e con la risposta all'impulso si modella
% l'equivalente FIR.


clear;      % Azzera il workspace di MATLAB
close all;
FS=8000;


Fc=400;     %Frequenza di taglio del filtro IIR da implementare come FIR
a=1-exp(-(2*pi*Fc/FS));
b=1-a;

% Calcolo della risposta all’impulso h(n) del filtro y(n)=ax(n)+by(n-1)

h(1)=a;
n=2;
while 1
    h(n)=b*h(n-1);
    if h(n)<0.000001 break
    end
    n=n+1;
end

figure;
plot(h);

% Il seguente codice consente di visualizzare la risposta in frequenza del filtro IIR

H=fft(h,512); %(FFT Fast Fourier Tranform è la versione veloce della DFT)
RH=real(H); % componente reale di H
IH=imag(H); % componente immaginaria di H
AMP=sqrt(RH.^2+IH.^2)/512; % spettro di ampiezza
f = FS*(0:256)/512; % calibrazione asse frequenziale [Hz]
plot(f,AMP(1:257))
title('Frequency content of h')
xlabel('frequency [Hz]')

TAPS=n;  % Durata della risposta all’impulso

% Preparazione del test pattern per la visualizzazione della risposta in frequenza del filtro FIR

D=0.1*FS;
f=50;
for n=1:TAPS
    x(n)=0;
end
k=n;
while f<FS/2
    for n=1:D
        x(k)=sin(2*pi*n*f/FS);
        k=k+1;
    end
    f=f*2;
end
figure;
plot(x);

% Implementazione del filtro FIR
N=size(x);
n=TAPS+1;
while n<N(2)
    y(n)=0;
    for k=1:TAPS
       y(n)=y(n)+h(k)*x(n-k);
    end
    n=n+1;
end
figure;
plot(y);    % Visualizza la risposta in frequenza del filtro FIR (da comparare con quella del filtro IIR)
