Laboratorul 9.

Convoluția, filtre FIR și metoda de proiectare folosind ferestre

În acest laborator vom face câteva exerciții pentru a ne familiariza cu operația de convoluție precum și cu filtre cu răspunsul finit la impuls (finite impulse response - FIR) și cu metode de proiectare a acestora prin metode folosind ferestre.

Materiale utile:

  • Vedeți slide-urile de la R. Lyons aici

Exercițiul 1 -- Convoluția [5p]

Precum am spus la curs, putem defini operația de convoluție dintre două secvențe h(k) și x(n) după cum urmează:

$y(n) = h(k) \ast x(n) = \sum_{k=0}^{M-1} h(k) \cdot x(n-k)$,

unde $M$ este lungimea secvenței h(k). De reținut că această operație definește un singur element de ieșire y(n). Pentru următorul element de ieșire, y(n+1), trebuie sa shiftăm secvența h(k) astfel încât să se potrivească cu elementele [x(n+1), x(n), …, x(n-M+2)]. Vedeți slide-urile 5 și 6 aici.

În general presupunem că secvența h(k) are un număr $M$ finit și, în general, mic de elemente.

Rezolvați următoarele exerciții:

  1. Dacă x(n) are N elemente și h(k) are M elemente, câte elemente are secvența obținută prin convoluție $h(k) \ast x(n)$? (presupunând că efectuăm convoluția doar pe elementele diferite de zero, adică atunci când elementele celei mai scurte secvențe se suprapun în totalitate cu elementele celeilalte)
  2. Fie x(n) secvența [1, 3, 5, 7, 5, 4, 2] și h(k) secvența [0.1, 0.3, 0.1]. Secvența obținută prin convoluție $y(n) = h(k) \ast x(n), n \in \{1,2,3,4,5\}$ ar trebui sa aibă 5 elemente. Scrieți fiecare dintre aceste elemente ca un produs scalar. Observați că trebuie să inversați ordinea secvenței x(n) (doar partea care se înmulțeste cu filtrul) înainte să o înmulțiți cu h(k) (alternativ, puteți inversa o singură dată elementele filtrului).
  3. Fie x(n) o secvența de $N=64$ elemente corespunzătoare unei sinusoide de frecvență $f=3$ kHz, eșantionată cu $f_s=64$ kHz. Fie h(k) secvența [0.1, 0.2, 0.2, 0.2, 0.1]. Generați aceste secvențe în Matlab și implementați convoluția pentru a obține elementele $y(n) = h(k) \ast x(n)$. Plotați inputul x(n) și ieșirea y(n); folosiți funcția stem în locul funcției plot pentru acest exercițiu.
  4. Înlocuiți x(n) cu secvența de impuls cu 9 elemente [0, 0, 0, 0, 1, 0, 0, 0, 0] și efectuați convoluția cu aceeași secvență h(k) ca mai sus. Ce obțineți ca y(n)? Cum se numește aceasta?
  5. Încercați operațiile de mai sus folosind funcția Matlab conv. Obțineți aceleași rezultate? Care sunt diferențele?

Pentru a plota secvențe de lungimi diferite în același plot, doar afișați primele N elemente ale secvenței, renunțând la restul.

O posibilă soluție:

convolution.m
clear;
close all;
 
%% Generate the sinewave signal    
fs = 64000;
N = 64;
x = 0:(N-1);
f = 3000;
A = 1;
xn = A*sin(2*pi*f/fs*x);
h = figure;
stem(x, xn);
xlabel('Sample index');
title('Input sinusoid sequence');
 
%% Generate the coefficients h(k)
hk = [0.1, 0.2, 0.2, 0.2, 0.1];
M = length(hk);
 
%% Convolve h with x
% Note: we must invert the sequence x,
% here done by indexing with -1
L = N-M+1;
y = zeros(1, L);
for i=1:L
    y(i) = hk * xn(i+M-1:-1:i)'; % vector product
end
y2 = conv(xn, hk);
h = figure;
stem(1:L, xn(1:L), 'k');
hold on;
stem(1:L, y, 'b');
stem(1:L, y2(1:L), 'r');
xlabel('Sample index');
legend('x(n)', 'convolution only over non-zero', 'Matlab conv');
title('y(n) = x(n) convolved with h(k)');

Exercițiul 2 -- Filtre FIR [5p]

În acest exercițiu vom creea secvența x(n) pentru un filtru trece-jos plecând de la un filtru ideal în domeniul frecvență și, folosind IFFT, vom obține secvența x(n). Vom folosi de asemenea diferite ferestre pentru a le compara performanța în proiectarea de filtre trece-jos.

Urmați următorii pași:

  1. Generați o secvență de filtru ideal trece-jos HK având N = 256 elemente, reprezentând spectrul de frecvență al unui filtru trece-jos. Folosiți o frecvență de cut-off de fs/16. Adică totul înainte de fs/16 trebuie sa treacă, pe când totul mai sus trebuie sa fie oprit (folosiți un dreptunghi care se opreste la fs/16). Observați că trebuie să generați un spectru simetric pentru a obține o secvență reală la următorul pas. Ar trebui să obțineți ceva precum: Țineți minte că acesta poate fi văzută ca ieșirea din DFT(FFT), adică primul element corespunde frecvenței 0, pe când următoarele N/2 -1 corespund frecvențelor pozitive, iar ultimele N/2 componente reprezintă frecvențele negative.
  2. Acum aplicați inversa DFT(în practică inversa FFT, ifft în Matlab) pentru a obține secvența corespunzătoare în domeniul timp hk(n). Rețineți: trebuie sa aplicați funcția ifftshift pe rezultat pentru a obține o funcție sinc simetrică, adică folosiți ceva precum:
    hk = ifftshift(ifft(HK));
  3. Trunchiați secvența hk(n) prin selectarea a doar L=65 de eșantioane din centru(32 din stânga maximului funcției sinc, maximul funcției, și 32 de eșantioane din dreapta). Aceasta corespunde multiplicării secvenței hk(n) cu o fereastră rectangulară centrată în punctul maxim al funcției sinc. Plotați secvența.
  4. Aplicați DFT(fft) pe secvența trunchiată și plotați spectrul. Rețineți: este important aici, precum și la primul plot pentru filtru trece-jos ideal, să notăm axa frecvențelor (axa x) ca o funcție de Fs, adică de la 0 la 1. Vedeți diferențe fața de un filtru ideal trece-jos. Acestea sunt efectele ferestrei dreptunghiulare.
  5. Folosiți aceeași secvență trunchiată ca mai sus dar înmulțiți-o cu o fereastră precum Blackman (blackman in MATLAB). Efectuați din nou DFT și plotați spectrul. Arată mai bine?.
  6. În final, folosiți ca intrare sinusoida din Exercițiul 1 ca x(n) și filtrați-o printr-o convoluție cu secvența obținută mai sus după folosirea ferestrei Blackman(folositi functia conv din Matlab). Plotați intrarea și ieșirea în aceeași figură folosind stem pentru a observa efectele filtrului.

O posibilă soluție:

lpf_window.m
clear;
close all;
 
%% Generate an ideal-filter sequence (in freq domain)
N = 256;
kc = N/16;
HK = zeros(1, N);
HK(1:kc) = 1;
HK(end-kc+2:end) = 1;
h = figure;
fx = linspace(0,1,N);
plot(fx, HK);
xlabel('Frequency (x Fs)');
title('Ideal low-pass filter');
print(h, '-dpng', 'ideal_filter.png');
 
%% Obtain time-domain sequence
hk = ifftshift(ifft(HK));
figure;
plot(1:N, hk);
xlabel('Index');
title('Time-domain sequence of ideal low-pass filter');
 
%% Generate windowed time-domain sequence
% window = rect
L = 65;
w = ones(1, L);
hw = hk(N/2+1-floor(L/2):N/2+1+floor(L/2)) .* w;
h = figure;
plot(1:L, hw);
xlabel('Sample index');
title('Windowed filter sequence h(k) using rect(n)');
 
%% Obtain frequency-domain of windowed sequence
HW = fft(hw);
figure;
fx = linspace(0,1,L);
plot(fx, abs(HW));
xlabel('Frequency (x Fs)');
title('FFT of the time-domain sequence with w(n)=rect(n)');
 
%% Generate windowed time-domain sequence
% window = blackman
L = 65;
w = blackman(L)';
hw = hk(N/2+1-floor(L/2):N/2+1+floor(L/2)) .* w;
h = figure;
plot(1:L, hw);
xlabel('Sample index');
title('Windowed filter sequence h(k) using blackman(n)');
 
%% Obtain frequency-domain of windowed sequence
HW = fft(hw);
figure;
fx = linspace(0,1,L);
plot(fx, abs(HW));
xlabel('Frequency (x Fs)');
title('FFT of the time-domain sequence with w(n)=blackman(n)');
 
%% Use the sequence to filter a sinewave
fs = 64000;
N = 64;
x = 0:(N-1);
f = 3000;
A = 1;
xn = A*sin(2*pi*f/fs*x);
yn = conv(xn, hw);
h = figure;
stem(x, xn, 'b');
hold on;
stem(x, yn(1:N), 'r');
xlabel('Sample index');
legend('xn', 'yn');
title('Filtering using hw');

ps/laboratoare/09.txt · Last modified: 2018/12/08 21:14 by ionut.gorgos
CC Attribution-Share Alike 3.0 Unported
www.chimeric.de Valid CSS Driven by DokuWiki do yourself a favour and use a real browser - get firefox!! Recent changes RSS feed Valid XHTML 1.0