CPS

konfrontacja_filtrow_IIR.m

% m-plik skrytpowy: rozwiazanie_zadania_3.m
%
% Laboratorium: Cyfrowe Przetwarzanie Sygnałów;
% Katedra Elektroniki, WIEiT, AGH;
%
% Opracowanie: Przemyslaw Koerohoda,, 12/05/2015;

clc; clear; close all;

fp=11025; fg=1000;
M=50;
h0=fir1(M,fg/(fp/2));
%h0 = rand(1,M+1);
h1=h0.*(-1).^(0:M); % fn = fp - fg
d=zeros(1,M+1); d(M/2+1)=1; % fn = fg
h2=d-h0;

    %wykresy_filtru_1(h0,1,1000,fp,1,M+1);
    %wykresy_filtru_1(h1,1,1000,fp,2,M+1);
    %wykresy_filtru_1(h2,1,1000,fp,3,M+1);
    wykresy_filtru_1(h0+h1,1,1000,fp,4,M+1);
    wykresy_filtru_1(d-(h0+h1),1,1000,fp,5,M+1);
% KONIEC PLIKU;

korohoda1.m

% m-plik skrytpowy: rozwiazanie_zadania_3.m
%
% Laboratorium: Cyfrowe Przetwarzanie Sygnałów;
% Katedra Elektroniki, WIEiT, AGH;
%
% Opracowanie: Przemyslaw Koerohoda,, 12/05/2015;

clc; clear; close all;

fp=11025; fg=1000;
M=50;
h0=fir1(M,fg/(fp/2));
%h0 = rand(1,M+1);
h1=h0.*(-1).^(0:M); % fn = fp - fg
d=zeros(1,M+1); d(M/2+1)=1; % fn = fg
h2=d-h0;

    %wykresy_filtru_1(h0,1,1000,fp,1,M+1);
    %wykresy_filtru_1(h1,1,1000,fp,2,M+1);
    %wykresy_filtru_1(h2,1,1000,fp,3,M+1);
    wykresy_filtru_1(h0+h1,1,1000,fp,4,M+1);
    wykresy_filtru_1(d-(h0+h1),1,1000,fp,5,M+1);
% KONIEC PLIKU;

krohoda2.m

% m-plik skryptowy: demo_filtru_IIR_niestabilnego.m
%
% Laboratorium: Cyfrowe Przetwarzanie Sygnałów;
% Katedra Elektroniki, WIEiT, AGH;
%
% Opracowanie: Przemyslaw Korohoda, 19/05/2015;


clc; clear; close all;
fp=1000;

opcja=3;
switch opcja,
    case 1, K=3; p=rand(K,1)+rand(K,1)*j; p=[p;conj(p)]; s=max(abs(p)); p=p/s*0.9;
    case 2, p=[1+j;1-j;0.5;-0.5];
    case 3, K=3; p=rand(K,1)+rand(K,1)*j; p=[p;conj(p)]; s=max(abs(p)); p=p/s*1.1; 
end;

M=length(p);
b=rand(1,M); a=poly(p);
b=b/a(1); a=a/a(1);

f1=33; N=2^10; dt=1/fp;
t=0:dt:(N-1)*dt;
x=cos(2*pi*f1*t);
y=filter(b,a,x);

[Hx,fx]=freqz(b,a,[f1,2*f1],fp),
MHx=abs(Hx),
PHx=angle(Hx),
y1=MHx(1)*cos(2*pi*f1*t+PHx(1));

    figure(1);
        freqz(b,a,1000,fp);
        subplot(2,1,1); hold on; plot(f1,20*log10(MHx(1)),'r.');
        subplot(2,1,2); hold on; plot(f1,PHx(1)/pi*180,'r.');
    figure(2);
        zplane(b,a);
    figure(3);
        subplot(2,1,1); plot(t,x,'b.-'); grid on;
                             xlim([t(1),t(end)]);
        subplot(2,1,2); plot(t,y,'g.-'); grid on; hold on;
                        plot(t,y1,'r.');
                             xlim([t(1),t(end)]); ylim([-2*MHx(1),2*MHx(1)]);

% KONIEC PLIKU

optym_IIR.m

function [b,a]=optym_IIR(M,fg,fp);
% [b,a]=optym_IIR(M,fg,fp);
%
% Laboratorium: Cyfrowe Przetwarzanie Sygnałów;
% Katedra Elektroniki, WIEiT, AGH;
%
% Opracowanie: Przemyslaw Korohoda, 19/05/2015;


N=1000;
f_=(0:(N-1))/(2*N);
MHref=zeros(1,N);
MHref(f_<=fg/fp)=1;
Href=MHref.*exp(-j*(2*pi*(f_)*M/2));
ww=ones(1,N);
%[b0,a0]=ellip(M,3,60,(fg/(fp/2)));
[b0,a0]=butter(M,(fg/(fp/2)));
ba0=[b0,a0];

      Opcje=optimset('TolFun',1e-14,'MaxFunEvals',100000,'MaxIter',100000);
ba=fminsearch(@szukam_filtru_IIR,ba0,Opcje,Href,ww);
b=ba(1:M+1);
a=ba(M+2:end);

% KONIEC FUNKCJI;q

szukam_filtru_IIR.m

function  F=szukam_filtru_IIR(ba,Href,ww);
% F=szukam_filtru_IIR(ba,Href,ww);
%
% Laboratorium: Cyfrowe Przetwarzanie Sygnałów;
% Katedra Elektroniki, WIEiT, AGH;
%
% Opracowanie: Przemyslaw Korohoda, 19/05/2015;

N=length(Href);
M=length(ba)/2-1;
b=ba(1:M+1);
a=ba(M+2:end);
[H,f]=freqz(b,a,N,1);

F=sum(abs( H(:)-Href(:) ).^2);

% KONIEC FUNKCJI; 

test_stabilnosci_wg_Hurwitza.m

% m-plik skryptowy: test_stabilnosci_wg_Hurwitza.m
%
% Laboratorium: Cyfrowe Przetwarzanie Sygnałów;
% Katedra Elektroniki, WIEiT, AGH;
%
% Opracowanie: Przemyslaw Korohoda, 19/05/2015;


clc; clear; close all;

opcja=3;
switch opcja,
    case 1, K=3; p=rand(K,1)+rand(K,1)*j; p=[p;conj(p)]; s=max(abs(p)); p=p/s*0.9;
    case 2, p=[1+j;1-j;0.5;-0.5];
    case 3, K=3; p=rand(K,1)+rand(K,1)*j; p=[p;conj(p)]; s=max(abs(p)); p=p/s*1.1; 
end;

M=length(p);
b=rand(1,M); a=poly(p);
b=b/a(1); a=a/a(1);

% przenosimy sie do dziedziny "s":
ps=(p-1)./(p+1);
as=poly(ps); as=as/as(1),

Hurwitz=zeros(M,M);
as_inv=[zeros(1,M+1),as(end:-1:1),zeros(1,M+1)],
for m=1:M,
  Hurwitz(m,:)=as_inv(((2*(M+1)-1):(3*(M+1)-3))-2*(m-1));  
end;

Hurwitz,

% wyznaczniki (minory wiodace):
for m=1:M,
   Det(m)=det(Hurwitz(1:m,1:m)); 
end;
Det,

    figure(1);
        subplot(1,2,1); plot(real(p),imag(p),'bo'); hold on;
                        plot(sin(0:0.001:3*pi),cos(0:0.001:3*pi),'k');
                         axis equal; grid on;
        subplot(1,2,2); plot(real(ps),imag(ps),'bo'); grid on;

% KONIEC PLIKU;

wykresy_filtru_1.m

function wykresy_filtru_1(b,a,N,fp,fn,K);
% wykresy_filtru_1(b,a,N,fp,fn,K);
%
% Funkcja rysuje wykresy opisujšce filtr z czasem dyskretnym;
% Nie zawiera opisu osi itd. - mozna te elementy uzupelnic dla poprawy czytelnosci wykresow.
%
% Laboratorium: Cyfrowe Przetwarzanie Sygnałów;
% Katedra Elektroniki, WIEiT, AGH;
%
% Opracowanie: Przemyslaw Korohoda, 5/05/2015;

[H,f]=freqz(b,a,N,fp);
h=impz(b,a,K);
Hr=exp(-j*(2*pi*(f/fp)*(K-1)/2));

    figure(fn); clf;
        subplot(2,2,1); plot(f,abs(H),'b.-'); grid on; xlim([0,fp/2]);
        subplot(2,2,2); plot(f,20*log10(abs(H)),'b.-'); grid on; xlim([0,fp/2]);
        subplot(2,2,3); plot(0:K-1,h,'k.-'); grid on; xlim([0,K]);
        subplot(2,2,4); plot(f,angle(H)/pi,'b.-'); grid on; xlim([0,fp/2]); hold on;
                        plot(f,angle(Hr)/pi,'r.');
                        
% KONIEC FUNKCJI;

wykresy_filtru_2.m

function wykresy_filtru_2(b,a,N,fp,fn,K,M);
% wykresy_filtru_2(b,a,N,fp,fn,K,M);
%
% Funkcja rysuje wykresy opisujšce filtr z czasem dyskretnym;
% Nie zawiera opisu osi itd. - mozna te elementy uzupelnic dla poprawy czytelnosci wykresow.
%
% Laboratorium: Cyfrowe Przetwarzanie Sygnałów;
% Katedra Elektroniki, WIEiT, AGH;
%
% Opracowanie: Przemyslaw Korohoda, 18/05/2015;


[H,f]=freqz(b,a,N,fp);
h=impz(b,a,K);
Hr=exp(-j*(2*pi*(f/fp)*M/2));

    figure(fn); clf;
        subplot(2,2,1); plot(f,abs(H),'b.-'); grid on; xlim([0,fp/2]);
        subplot(2,2,2); plot(f,20*log10(abs(H)),'b.-'); grid on; xlim([0,fp/2]);
        subplot(2,2,3); plot(0:K-1,h,'k.-'); grid on; xlim([0,K]);
        subplot(2,2,4); plot(f,angle(H)/pi,'b.-'); grid on; xlim([0,fp/2]); hold on;
                        plot(f,angle(Hr)/pi,'r.');
                        

% KONIEC FUNKCJI;