FIR DP

%skrypt TypI_DP.m
%projektowanie filtru FIR DP o liniowej fazie
%metoda probkowanie w dziedzinie czestotliwosci (i okna)
%Typ 1 - symetryczna odp,. impulsowa; liczba wspolczynnikow N=2M+1
%Typ 1 - czesc urojona funkcji transmitancji IM[H(m)]=0
% w tym przykładzie czestotliwosc graniczna filtru DP, fgr = 0.25*fpr
%%
clear all;clc; clf; subplot(1,1,1)
%N=33;
%N=65;
N=129; %RZAD FILTRU
disp(['Rzad filtru N=', num2str(N)])

% tworzenie idealnej f. transmitancji (1-pasmo przepustowe; 0 -p asemo
% zaporowe)

M=(N-1)/2; % 0.5 rzedu filtru
M2=M/2; % lub to samo: M2=(Nn-1)/4
j2=ones(1,M2);
z2=zeros(1,M2);
Ar=[1 j2 z2 z2 j2]; % czesc rzeczywista funkcji transmitancji
Ai=zeros(1,N); %czesc urojona funkcji transmitacji (zeruje sygnał)
A=Ar+1j*Ai; % transmitancja idealna (zespolone wartości)
%% - obliczanie IDFT (A) ----> h(n) <--odpowiedz impusowa
n=0:(N-1); % wektor indeksowy odpowiedzi impulsowej h(n)
f=n/N; % wektor wartosci czestotliwosci unormowanej
h=zeros(1,N); % baza sumowania
for k=0:(N-1)
    h=h+A(k+1)*exp(1j*2*pi*k/N*(n-M));
end
h = real(h)/N; % normowanie IDFT
okno=hamming(length(h));
h=h.*okno';
subplot(3,2,1)
stem(n-M,h); grid on
hold on
plot(n-(N-1)/2,okno'*max(h),'r')
title(['Odp. impulsowa h(n) filtru FIR. Rzad filtru N=',num2str(N)]);

%% =====Rzeczywista funkcja transmitancji H(m) ==============
NF=1000; % liczba probek transmitancji filtru
k=0:NF-1; % wektor indeksowy dla kolejnych wartosci H(m)
fn=k/NF; % wektor wartości czestotliwosci unormowanej
for k=0:NF-1 % DFT z definicji
    H(k+1)=exp(-1i*2*pi*k/NF*(n-N))*h';
    % to jest iloczyn skalarny wektora wartosci exp() i wektora h'
end
Ax=Ar; %funkcja transmitancji idealna
Hx=real(H); % funkcja transmitancji rzeczywista
subplot(3,2,3)
stem(f,Ax,'or'); grid on
hold on
plot(fn,abs(H),'b')
title('Charakterystyka filtru');
legend('idealna','faktyczna')
subplot(3,2,5)
plot(fn,20*log10(abs(H)));grid on
title('|H(m)| w [db]')
%%
%%=====tworzenie sygnalu testowego i jego charakterystyki
fs=1000;  %czestotliwosc probkowania [Hz]
deltaT=1/fs; %okres probkowania
t=0:deltaT:1-deltaT; %wektor chwili czasu probkowania (od 0 do 1 sekundy)
f1=100; f2=350; % czestotliwosc sinusoid
s1=sin(2*pi*f1*t); % 1. sinusoida
s2=sin(2*pi*f2*t); % 2. sinusoida
syg=s1+s2; % sygnal tesotwy (2 sinsoidy)
syg=syg+2*randn(1,length(syg)); % sygnal testowy z szumem
%%
% ykres
subplot(3,2,2)
plot(t,syg); grid on; xlabel('czas[s]')
title('Sygnal testowy')
F=fft(syg); %dft Syg
Fsyg=abs(F); %moduł z syg czyli widmo amplitudowe
subplot(3,2,4)
plot(Fsyg); grid on; title('Widmo amplitudowe sygnału testowego')
%%
% splot
% czyli filtracja sygnału
y=conv(syg,h);
FF=abs(fft(y));
subplot(3,2,6)
plot(FF); grid on; title('Widmo aplitudowe po filtracji')