Signalerx6


ØVELSE 6
i 31610


2009-04-06

Rasmus Berg Palm s062241


Mads Andersen s062237


Jacob Bjerring Olesen s062248


Indholdsfortegnelse

Excercise 1

Med MATLABs konvention er et lineært filter beskrevet ved vektorerne $A=[a_0, a_1, a_2, ..., a_N], B=[b_0, b_1, b_2, ..., b_M]$, et system der kan beskrives ved følgende lineære differens ligning med konstante koefficienter:

(1)
\begin{align} a_0\cdot y(n) = \sum^{M}_{k=0}b_k\cdot x(n-k) -\sum^{N}_{k=1}a_k\cdot y(n-k) \end{align}

Det bemærkes at for et sådan system, afhænger output ikke af fremtidige input, og systemet er derfor kausalt.
Systemet beskrives i z-domænet som:

(2)
\begin{aligned} a_0\cdot Y(z) &= \sum^{M}_{k=0}b_k\cdot X(z)z^{-k} -\sum^{N}_{k=1}a_k\cdot Y(z)z^{-k} \quad \Leftrightarrow\\ H(z) &=\frac{Y(z)}{X(z)}=\frac{\sum^{M}_{k=0}b_k\cdot z^{-k}}{\sum^{N}_{k=0}a_k\cdot z^{-k}} \end{aligned}

Sættes $z^{-M}$ uden for parentes i tælleren, og $z^{-N}$ uden for parentes i nævneren fås:

(3)
\begin{align} H(z) =\frac{z^{-M}}{z^{-N}} \frac{\sum^{M}_{k=0}b_{M-k}\cdot z^{k}}{\sum^{N}_{k=0}a_{N-k}\cdot z^{k}}=z^{-M+N}\frac{b_0\cdot z^M+b_1\cdot z^{M-1}+\cdots +b_M}{a_0\cdot z^N+a_1\cdot z^{N-1} + \cdots + a_N} \end{align}

Rødderne til polynomiet i tælleren, z=z_k, k=1,..,M er nulpunkter for systemet, rødderne til polynomiet i nævneren, z=p_k, k=1,..,N er poler for systemet. Derudover er der -M+N poler i z=0, hvis M>N, -M+N nuller i z=0 og en pol i $z=\inf$ , hvis N>M.

ROC for et stabilt LTI-system skal inkludere enhedscirklen (se bevis s196-197 i Proakis bog). ROC for et kausalt IIR-system er det ydre af en cirkel med radius r, dvs for et stabilt kausalt IIR-LTI-system skal r<1. For z=p_k divergerer H(z), dvs. et kausalt LTI-system er stabilt hvis og kun hvis |p_k|<1 (s197 i Proakis bog).

Funktionen isstable.m tjekker stabilitet af et filter givet ved vektorerne A og B der beskriver et lineært filter efter MATLABs konvention.
Roots.m bruges til at finde poler og nulpunkter for systemet, idet roots(A) finder rødderne til polynomiet $a_0\cdot z^N+a_1\cdot z^{N-1} + \cdots + a_N$. Der tjekkes derefter om en pol og et nulpunkt ophæver hinanden, og de resterende polers belligenhed tjekkes om er indefor enhedscirklen. Hvis polerne ligger indefor enhedscirklen er systemet stabilt og isstable returnerer stable=1.

function [stable,p,z] = isstable(a,b)
% Finds stability of a system described by the constant coefficient 
% difference eqaution:
% a(1)*y(n) = b(1)*x(n) + b(2)*x(n-1) + ... + b(nb+1)*x(n-nb)
%                       - a(2)*y(n-1) - ... - a(na+1)*y(n-na)
% If system is stable, istable.m gives stable=1

stable = 0;
p = roots(a); %the poles of H(z)
z = roots(b); %the zeroes of H(z)

%Pole-zero cancellation:
% if there is a pole same place as a zero, the pole and zero cancel out
% each other
for i=p 
    x=find(i==z); 
    p(x)=[];
    z(x)=[];
end

    %If the poles lay inside the unit circle, the system is stable:
if(abs(p)<=1) 
    stable = 1;
end

isstable blev brugt på filtrene 1, 2, 3 med resultat:

stable1 =

     1

stable2 =

     0

stable3 =

     0

filter 1 er altså stabilt, filter 2 og 3 er ustabile. Dette kan også ses ved at påvirke systemerne med en impuls og iagtage impulsresponsene. En impuls er derfor dannet som en vektor med bestående af et 1-tal efterfulgt af 64 0'er. Respons fra systemerne opnås ved brug af MATLAB funktionen filter.m der udregner respons fra systemerne ved brug af en "Direct Form II Transposed"-implementering. filter.m er derfor brugt på filtrene med impulsen som input, og impulsresponsene ses plottet nedenfor:

h1.png
h2.png
h3.png

Det ses at impulsrespons fra filter 2 og 3 divergerer, hvorimod impulsrespons fra filter 1 konvergere mod 0, hvilket klart bekræfter at filter 1 er stabilt og filter 2 og 3 er ustabile.

Følgende MATLAB kode blev benyttet i denne opgave:

clear all; close all; clc;
%% Ex 1
% systems are defined:
B1=[0.0725 0.2200 0.4085 0.4883 0.4085 0.2200 0.0725];
A1=[1.0000 -0.5835 1.7021 -0.8477 0.8401 -0.2823 0.0924];
B2=[1.0000 1.3000 0.4900 -0.0130 -0.0290];
A2=[1.0000 -0.4326 -1.6656 0.1253 0.2877];
B3=[1.0000 -1.4000 0.2400 0.3340 -0.1305];
A3=[1.0000 0.5913 -0.6436 0.3803 -1.0091];

%stability is determined using isstable.m:
[stable1, p1, z1] = isstable(A1,B1);
[stable2, p2, z2] = isstable(A2,B2);
[stable3, p2, z2] = isstable(A3,B3);

%impulse is defined:
delta = [1 zeros(1,64)];

%impulse responses are determined using filter.m:
h1 = filter(B1,A1,delta);
h2 = filter(B2,A2,delta);
h3 = filter(B3,A3,delta);

%plots of impulse responses:
figure;
plot(h1)
xlabel('n')
ylabel('h_1(n)')
title('Impulse response of filter 1')

figure;
plot(h2)
xlabel('n')
ylabel('h_2(n)')
title('Impulse response of filter 2')

figure;
plot(h3)
xlabel('n')
ylabel('h_3(n)')
title('Impulse response of filter 3')

Excercise 2

Vi har fra Excercise 1 at

(4)
\begin{align} H(z) &=\frac{Y(z)}{X(z)}=\frac{\sum^{M}_{k=0}b_k\cdot z^{-k}}{\sum^{N}_{k=0}a_k\cdot z^{-k}} \end{align}

Sættes $z = e^{j 2 \pi f}$ fås:

(5)
\begin{align} H(f) =\frac{\sum^{M}_{k=0}b_k\cdot e^{-j 2 \pi f k}}{\sum^{N}_{k=0}a_k\cdot e^{-j 2 \pi f k}} = \frac{B(f)}{A(f)} \end{align}

Det ses altså at tælleren er fouriertransformationen af vektoren B og at nævneren er den fouriertransformerede af A.
Dette kan implementeres effektivt i matlab vha. fft:

function [H, f] = transfer(A, B, nf)
% A, B are coefficient values defined as in filter.m. nf are number of
% points in the spectrum of calculation.
f = -0.5:1/(nf-1):0.5;

Y = fftshift(fft(B,nf));
X = fftshift(fft(A,nf));

H = Y./X;

Excercise 3

We apologize for the mixed-up language..

The stable filter given by B1 and A1 filtered a delta function of length 64. Afterwards it got transformed into the frequency domain by using the MATLAB function fft. In this domain the delta function becomes a transfer function, which is compared graphically with the one found in a previous question. At first the two graphs lookalike, except for the more fluctuating curve seen for the new transfer function.

trans.png
fourier.png

But in a more thorough view you notice a slightly resembling between the lobes and a sinc function. This might have something to do with the zero-padding taken place in the fft function, which is equivalent to multiplying the signal by a rectangular window. This multiplication is equivalent to convolving the transformed signal with a sinc function in the frequency domain.

%ex2:
[H1, f1] = transfer(A1,B1,100000);
figure;
plot(f1,db(abs(H1)))
ylim([-120 20])
title('The transfer function found by "transfer.m"')

%ex3

delta = [1 zeros(1,63)];%A delta function of length 64
df = filter(B1,A1,delta);%Filtering delta by the stable filter given in ex1
DF = fftshift(fft(df,100000));%Shifting df into the frequency domain
figure
plot(f1,db(abs(DF)))
ylabel('Magnitude dB')
xlabel('Normalized frequency')
title('The Fourier transformed of the impuls respons')

Excercise 4

En orden på 64 viste sig at være nødvendig for at opnå en 45 dB dæmpning ved 570 Hz. Nedenfor ses et plot af overføringsfunktionen

firpm.png

Følgende kode blev benyttet til denne opgave

pb = 500; %Passband frequency
sb = 570; %Stopband frequency
fs = 2000; %Sampling frequency
fnq = fs/2; %Nuqvuist frequency
n = 64; %Filter order

f = [0 pb/fnq sb/fnq 1]; %Important filter Frequencies
a = [1 1 0 0]; %Amplitude corresponding to values of f
b = firpm(n,f,a); %Creates the filter coefficients
[H, w] = freqz(b,1,512); %Works similar to transfer
sb_att = db(H(find(w/pi*fnq>570,1,'first'))); %Finds the stopband attenuation
figure
hold on;
plot(w/pi*fnq,db(abs(H))) %Plot the filter
h=line([500 500],[-100 20]); set(h,'Color','g'); %Plot passband line
h=line([570 570],[-100 20]); set(h,'Color','r'); %Plot stopband line
h=line([0 fnq],[sb_att sb_att]); set(h,'Color','k'); %Plot attenuation line
hold off
title('firpm transfer function');
ylabel('Attenuation [db]')
xlabel('Frequency [Hz]')
legend([ num2str(n) ' order firpm filter'],'Passband (500 Hz)','Stopband (570 Hz) ', [num2str(sb_att) ' dB attenuation'])

Excercise 5

In the class of IIR filters, elliptic filters are the most efficient to use in the sense that for some given specification the elliptic filter has the lowest order in fulfilling these, and thereby the fewest number of calculations. As the order increased the approximated error decreased, which makes the filter more precise on the count of more calculations. The order of a digital filter is a number describing the highest exponent in the numerator or denominator of the z-domain transfer function of a digital filter. For IIR filters, the filter order is equal to the number of delay elements in the filter structure.
As seen by our routine, a proper filter order is in the vicinity of 6 and then become more and more precise up to the order of 16.

ellip6
epilic16.png

In comparison to the previous transfer function, you’ll see that the order of the elliptic filter has decreased by a factor 10! This decrease results in fewer calculation, which is preferable in most software.
We have had a hard time estimating the specific number of calculation made by the filter. This is due to the fact that the amount of calculation done by a filter isn’t always the same, even though the order is. We therefore need to have more knowledge about the mathematical appearance of the transfer function, to suggest the number of calculation needed.
One of the only reason not to those an elliptic filter is due to its more nonlinearly phase response, compared to the response given by either a Butterworth- or Chebyshev filter

[[code]]
%ex5
%Designing a Lowpass filter with the ellip function.
%ellip(N,Rp,Rs,Wp)
%N = order of the filter
%Rp = The peak to peak ripple given in units of decibel (0.5dB)
%Rs = The minimum stop-band attenuation given in decibel (45dB)
%Wp = The passband-egde frequency. Wp must be 0.0 < Wp < 1.0,
% with 1.0 corresponding to half the sample rate (Wp = 0.5)
[B4,A4] = ellip(6,0.5,45,0.5,'low')

%Finding the coefficient of the transfer function
[H4, f4] = transfer(A4,B4,100000);
figure
plot(f4,db(abs(H4)))
ylim([-120 20])
ylabel('Magnitude dB')
xlabel('Normalized frequency')
title('Elliptic filter of order 6')

[B5,A5] = ellip(16,0.5,45,0.5,'low')

%Finding the coefficient of the transfer function
[H5, f5] = transfer(A5,B5,100000);
figure
plot(f5,db(abs(H5)))
ylim([-120 20])
ylabel('Magnitude dB')
xlabel('Normalized frequency')
title('Elliptic filter of order 16')
[[\code]]

Excercise 6

Med MATLAB lowpass filter design demoen er dannet lowpass-filtre med specifikationerne:
pass-band indtil 500Hz
Stop-band fra 600Hz
amplitude af ripple i pass-band: 3dB
amplitude af ripple i stop-band: 50dB

Og Matlab har derudfra selv beregnet den nødvendige orden af filteret. Der er dannet FIR-filtrene FIRPM, FIRLS og KAISER, og IIR-filtrene BUTTER, CHEBY1, CHEBY2, og ELLIP. Resultaterne ses nedenfor:

FIR-filtrene:

firpm1.png
firls.png
kaiser.png

IIR-filtrene:

butter.png
cheby1.png
cheby2.png
ellip.png

Det ses at IIR-filtrene opnås de ønskede specifikationer ved en lavere orden, især ekstremerne er ellipfiltret i forhold til kaiser filtret. Men igen har vi haft svært ved at bestemme den konkrete antal udregninger som filteret har gennemført, da vi som forklaret i forrige spørgsmål ikke kender overførings funktionens matematiske udseende.

edit sections

Beskeder
Kritik modtages på mine afsnit. VH, Ralle
kritik på ex1 modtages, vh Mads
"For z=p_k divergerer H(z)". p_k er ikke forklaret. vh Ralle
hva med lidt kommentarer til script i 4'eren
ikke minus, det er z=e^{j 2 \pi f}. Efter H(f) er opskrevet så måske lige sætningen: Det ses altså at tælleren er fouriertransformationen af vektoren B, nævneren er den fouriertransformerede af A...
sigx6
Here is a place for your title Click me to edit ! false
Unless otherwise stated, the content of this page is licensed under Creative Commons Attribution-ShareAlike 3.0 License