Podstawy analizy częstotliwościowej i próbkowanie sygnałów
Autorzy: Przemysław Korohoda, Krzysztof Duda
Pliki: rozdzial_02.zip
Przykład 2.1
Jest to program napisany w języku pakietu MATLAB, służący do wykreślania sygnału kosinusoidalnego (plik przyklad_2_1.m).% Przykład 2.1:
% m-plik: przyklad_2_1.m
%
% Cyfrowe przetwarzanie sygnałów. Podstawy, multimedia, transmisja. PWN, Warszawa, 2014.
% Rozdział 2. Podstawy analizy częstotliwościowej i próbkowanie sygnałów.
% Autorzy rozdziału: Przemysław Korohoda, Krzysztof Duda.
%
% Więcej na stronie internetowej: http://teledsp.kt.agh.edu.pl
%
% Opracowanie przykładu: kwiecień 2013 r./PK/KD.
% 1) przygotowanie środowiska MATLAB do uruchomienia programu:
clc; % czyścimy okno poleceń i przesuwamy kursor do lewego-górnego narożnika tego okna;
clear; % opróżniamy pamięć zmiennych;
close all; % zamykamy wszystkie okna graficzne;
% 2) przygotowanie najważniejszych parametrów:
dt=0.0001; % odstęp między wartościami na osi czasu, [s];
tmax=0.3; % maksymalna wartość dla przedziału czasu, [s];
f0=10; % wartość częstotliwości, [Hz];
A=3; % amplituda sygnału;
phi=pi/2; % faza początkowa [rad];
% 3) obliczenia:
t=0:dt:tmax; % ciąg kolejnych wartości na osi czasu, [s];
faza=2*pi*f0*t - phi; % wyznaczamy ciąg "faza", stała "pi" (=3,14...) jest wbudowana;
x=A*cos(faza); % korzystamy z funkcji "cos" - liczymy kolejno dla każdej wartości z ciągu "faza" wartości ciągu "x";
% 4) prezentacja graficzna:
figure(1); % otwieramy okno graficzne numer 1;
plot(t,x,'r.-'); % rysujemy wykres x(t) w kolorze czerwonym (r=red);
% kropki pokazują położenie punktów, z których powstał wykres "ciągły";
grid on; % na wykresie dodajemy siatkę linii poziomych i pionowych;
xlabel('t [s]'); % dodajemy opis osi poziomej;
% KONIEC PRZYKŁADU 2.1.
Przykład 2.2
W przykładzie (plik przyklad_2_2.m) pokazano jak w jednym oknie wytworzyć trzy wykresy (funkcja: subplot), na każdym wykresie umieścić sygnał, a następnie wrysować jego próbki (komenda: hold on). Ponadto opisujemy osie na wykresie (funkcje: xlabel, ylabel) i ograniczamy ich zakres w poziomie (funkcja: xlim).% Przykład 2.2:
% m-plik: przyklad_2_2.m
%
% Cyfrowe przetwarzanie sygnałów. Podstawy, multimedia, transmisja. PWN, Warszawa, 2014.
% Rozdział 2. Podstawy analizy częstotliwościowej i próbkowanie sygnałów.
% Autorzy rozdziału: Przemysław Korohoda, Krzysztof Duda.
%
% Więcej na stronie internetowej: http://teledsp.kt.agh.edu.pl
%
% Opracowanie przykładu: kwiecień 2013 r./PK.
clc; clear; close all;
dt=0.001; tmax=0.3; t=0:dt:tmax; % pseudo-ciągła oś czasu;
f1=10; f2=60; f3=12.1; % zadane częstotliwości w [Hz];
p1=pi/2; p2=pi/2; p3=pi/4; % zadane fazy początkowe w [rad];
A1=3; A2=3; A3=5; % % zadane amplitudy;
x1=A1*cos(2*pi*f1*t + p1); % kolejne linie są analogiczne - generują sygnały "ciągłe";
x2=A2*cos(2*pi*f2*t + p2);
x3=A3*cos(2*pi*f3*t + p3);
Dt=0.02; tn=0:Dt:tmax; % wybieramy punkty próbkowania dla x[n];
x1n=A1*cos(2*pi*f1*tn + p1); % trzy kolejne linie generują sygnały dyskretne;
x2n=A2*cos(2*pi*f2*tn + p2);
x3n=A3*cos(2*pi*f3*tn + p3);
figure(1); % część graficzna;
subplot(3,1,1); plot(t,x1,'r'); grid on; hold on;
plot(tn,x1n,'bo'); ylabel('x1'); xlim([0,tmax]);
subplot(3,1,2); plot(t,x2,'r'); grid on; hold on;
plot(tn,x2n,'bo'); ylabel('x2'); xlim([0,tmax]);
subplot(3,1,3); plot(t,x3,'r'); grid on; hold on;
plot(tn,x3n,'bo'); xlabel('t [s]'); ylabel('x3');
xlim([0,tmax]);
% KONIEC PRZYKŁADU 2.2.
Przykład 2.3
Przykład (plik przyklad_2_3.m)pokazuje, jak otrzymać graficzne oszacowanie ciągłej transformaty Fouriera zadanego sygnału o ograniczonym czasie trwania. Ponadto pokazuje: a) definiowanie własnej funkcji oraz korzystanie z tej funkcji, b) wykorzystanie funkcji szczególnie przydatnych dla liczb zespolonych (exp, abs, angle, real, imag), c) korzystanie z funkcji sinc, która zawiera wewnątrz mnożenie podanego argumentu przez stałą, co powoduje, że zapis w programie rożni się od zapisu według wzoru (2.6).% Przykład 2.3:
% m-plik: przyklad_2_3.m
%
% Cyfrowe przetwarzanie sygnałów. Podstawy, multimedia, transmisja. PWN, Warszawa, 2014.
% Rozdział 2. Podstawy analizy częstotliwościowej i próbkowanie sygnałów.
% Autorzy rozdziału: Przemysław Korohoda, Krzysztof Duda.
%
% Więcej na stronie internetowej: http://teledsp.kt.agh.edu.pl
%
% Opracowanie przykładu: kwiecień 2013 r./PK.
clc; clear; close all;
dt=0.001; tmax=300*dt; t=0:dt:tmax; % pseudo-ciągła oś czasu;
Nt=length(t); x=zeros(Nt,1); % przygotowujemy sygnał zapełniony wartościami zero;
T=tmax/2; N=T/dt-1;
x0=triang(N); x(2:N+1)=x0; % wstawiamy do sygnału trójkątny fragment;
a=T/2; t0=a; % parametry widma na podstawie teorii;
fmax=30; df=0.01; f=-fmax:df:fmax;
Nf=length(f); kk=1:Nf;
X1f(kk)=a*sinc(f*a).^2.*exp(-j*2*pi*f(kk)*t0); % widmo wyliczone na podstawie teorii;
for k=kk,
X2f(k)=prosta_calka(x'.*exp(-j*2*pi*f(k)*t),t); % widmo wyliczone z danych liczbowych;
end;
figure(1);
plot(t,x,'b.-'); grid on; xlabel('t [s]'); title('x(t)');
xlim([min(t),max(t)]);
figure(2);
subplot(2,2,1);
plot(f,abs(X1f),'r.-'); grid on; xlabel('f [Hz]'); axis tight; title('abs(X1(f))');
subplot(2,2,3);
plot(f,angle(X1f)/pi,'c.-'); grid on; xlabel('f [Hz]'); axis tight; title('angle(X1(f))/pi');
subplot(2,2,2);
plot(f,real(X1f),'g.-'); grid on; xlabel('f [Hz]'); axis tight; title('real(X1(f))');
subplot(2,2,4);
plot(f,imag(X1f),'m.-'); grid on; xlabel('f [Hz]'); axis tight; title('imag(X1(f))');
figure(3);
subplot(2,2,1);
plot(f,abs(X2f),'r.-'); grid on; xlabel('f [Hz]'); axis tight; title('abs(X2(f))');
subplot(2,2,3);
plot(f,angle(X2f)/pi,'c.-'); grid on; xlabel('f [Hz]'); axis tight; title('angle(X2(f))/pi');
subplot(2,2,2);
plot(f,real(X2f),'g.-'); grid on; xlabel('f [Hz]'); axis tight; title('real(X2(f))');
subplot(2,2,4);
plot(f,imag(X2f),'m.-'); grid on; xlabel('f [Hz]'); axis tight; title('imag(X2(f))');
% KONIEC PRZYKŁADU 2.3.
Przykład 2.4
Przykład (plik przyklad_2_4.m) pokazuje, jak na podstawie wzoru (2.14) wyznaczyć i wykreślić transformatę DTFT ciągu wyznaczonego w wyniku próbkowania sygnału z przykładu 2.3.% Przykład 2.4:
% m-plik: przyklad_2_4.m
%
% Cyfrowe przetwarzanie sygnałów. Podstawy, multimedia, transmisja. PWN, Warszawa, 2014.
% Rozdział 2. Podstawy analizy częstotliwościowej i próbkowanie sygnałów.
% Autorzy rozdziału: Przemysław Korohoda, Krzysztof Duda.
%
% Więcej na stronie internetowej: http://teledsp.kt.agh.edu.pl
%
% Opracowanie przykładu: kwiecień 2013 r./PK.
clc; clear; close all;
% najpierw tworzymy sygnał "analogowy" xa:
dt=0.001; tmax=300*dt; t=0:dt:tmax; Nt=length(t); xa=zeros(Nt,1);
N1=tmax/(2*dt); x0=triang(N1); xa(1:N1)=x0;
% następnie próbkujemy:
K=20; Dt=K*dt; % wyznaczamy odstęp próbkowania jako całkowitą wielokrotność dt (bo tak jest najwygodniej);
x=xa(1:K:end); % próbkujemy;
N=length(x); % wyznaczamy długość ciągu próbek;
% i poddajemy analizie Fouriera:
dF=0.0001; F=-0.5:dF:0.5; NF=length(F); % określamy pseudo-ciągły odcinek na osi F;
X=zeros(1,NF); % inicjalizacja pętli - wstępne wyzerowanie DTFT;
for n=0:N-1,
X=X+x(n+1)*exp(-j*2*pi*F*n); % dodajemy do DTFT części pochodzące od kolejnych x[n];
end;
font_size_1=14; % zmiana tej wartości zmieni wielkość czcionki w wybranych elementach opisu;
figure(1); % tę ilustrację warto rozciągnąć na cały ekran;
subplot(3,1,1);
plot((N-1)*t/tmax, xa, 'r'); grid on; hold on;
stem(0:N-1,x);
title('x[n]','fontsize',font_size_1); xlabel('n','fontsize',font_size_1);
subplot(3,2,3);
plot(F,abs(X),'r.-'); grid on;
xlim([-0.5,0.5]); title('abs(X)','fontsize',font_size_1);
subplot(3,2,5);
plot(F,angle(X)/pi,'c.-'); grid on;
title('angle(X)/pi','fontsize',font_size_1); xlabel('F','fontsize',font_size_1);
xlim([-0.5,0.5]);
subplot(3,2,4);
plot(F,real(X),'g.-'); grid on;
title('real(X)','fontsize',font_size_1);
xlim([-0.5,0.5]);
subplot(3,2,6);
plot(F,imag(X)/pi,'m.-'); grid on;
title('imag(X)','fontsize',font_size_1); xlabel('F','fontsize',font_size_1);
xlim([-0.5,0.5]);
% KONIEC PRZYKŁADU 2.4.
Przykład 2.5
W przykładzie 2.5 (plik przyklad_2_5.m) pokazano jak: a) zorganizować eksperyment komputerowy, pokazujący efekt aliasingu, b) zapisać sygnał akustyczny do pliku w formacie wave, c) odtworzyć na komputerze dźwięk zapisany w sygnale, d) sporządzić wykres dla wybranego fragmentu sygnału (wyświetlenie całości byłoby raczej nieefektywne – można spróbować!).% Przykład 2.5:
% m-plik: przyklad_2_5.m
%
% Cyfrowe przetwarzanie sygnałów. Podstawy, multimedia, transmisja. PWN, Warszawa, 2014.
% Rozdział 2. Podstawy analizy częstotliwościowej i próbkowanie sygnałów.
% Autorzy rozdziału: Przemysław Korohoda, Krzysztof Duda.
%
% Więcej na stronie internetowej: http://teledsp.kt.agh.edu.pl
%
% Opracowanie przykładu: kwiecień 2013 r./PK.
clc; clear; close all;
fpr=22025; dt=1/fpr; tmax=0.5; t=0:dt:tmax; % próbkowana oś czasu;
f1=440, f2=2*440, % warto dla porównania podstawić: f2=fpr-f1,
x1=cos(2*pi*f1*t);
x2=cos(2*pi*f2*t);
n0=1000; n=(1:100)+n0; % wybór indeksów na potrzeby wykresu;
figure(1); % początek części graficznej;
plot(t(n),x1(n),'bs-'); grid on; hold on;
plot(t(n),x2(n),'ro-');
xlabel('t [s]'); % opis osi poziomej;
axis tight; % ograniczamy wykres do "istotnego" zakresu;
x=[x1,x2,x1,x2,x1,x2]; % łączymy kolejno elementy sygnału;
wavwrite(x,fpr,16,'moja_muzyka_1.wav'); % zapis sygnału do pliku (wczytanie przez "wavread");
sound(x,fpr); % przesyłamy sygnał do karty dźwiękowej;
% KONIEC PRZYKŁADU 2.5.
Przykład 2.6
W przykładzie (plik przyklad_2_6.m) pokazano, jak zweryfikować oraz zilustrować dolnopasmową wersję twierdzenia o próbkowaniu dla zadanego sygnału x(t). Korzystamy z wzajemnych zależności między CFT oraz DTFT, które są opisane pod koniec rozdziału.% Przykład 2.6:
% m-plik: przyklad_2_6.m
%
% Cyfrowe przetwarzanie sygnałów. Podstawy, multimedia, transmisja. PWN, Warszawa, 2014.
% Rozdział 2. Podstawy analizy częstotliwościowej i próbkowanie sygnałów.
% Autorzy rozdziału: Przemysław Korohoda, Krzysztof Duda.
%
% Więcej na stronie internetowej: http://teledsp.kt.agh.edu.pl
%
% Opracowanie przykładu: kwiecień 2013 r./PK.
clc; clear; close all;
dt=0.0001; tmax=0.5; t=-tmax:dt:tmax; % odcinek pseudo-ciągłej osi czasu;
t=t(:); % ułożenie elementów ciągu t w wektorze kolumnowym;
Dt=0.0025; fpr=1/Dt; % Dt to faktyczny odstęp próbkowania;
tn0=0:Dt:tmax; tn=[-tn0(end:-1:2),tn0];% ciąg wartości t dla próbek, tak żeby w środku było t=0;
dF=0.0001; F=-2:dF:2; F=F(:); % odcinek pseudo-ciągłej osi F, F(:) - pionowo;
f=F*fpr; % odcinek pseudo-ciągłej osi częstotliwości f, w Hz;
f0=100; t0=0.01; % parametry sygnału;
x=sinc(f0*(t-t0)).^2; % sygnał pseudo-ciągły: x(t);
xn=sinc(f0*(tn-t0)).^2; % ciąg próbek: x[n];
NF=length(F); H=zeros(NF,1); % wstępne przygotowanie filtru odtwarzającego;
H((f>-fpr/2)&(f>fpr/2))=1; % dla wybranego zakresu f ustawiamy: H=1;
N=length(xn); X=zeros(NF,1); % przygotowanie do pętli liczącej DTFT;
for n=1:N,
n1=(n-(N+1)/2); % przeliczenie indeksów, by pasowały do wzoru DTFT;
X=X+xn(n)*exp(-j*2*pi*F*n1);
end;
X0=X.*H; % odtwarzamy widmo sygnału;
Nt=length(t); x_rec=zeros(Nt,1); % przygotowanie do pętli odtwarzającej sygnał;
fg=0.5*fpr; % częstotliwość graniczna filtru odtwarzającego;
for n=1:N,
n1=(n-(N+1)/2); % przeliczenie indeksów, by pasowały do wzoru;
t1=t-n1*Dt; % przesunięcie czasu dla każdego elementu sumy;
x_rec=x_rec+xn(n)*(2*fg/fpr)*sinc(2*fg*t1);
end;
rec_err=max(abs(x_rec-x)), % weryfikacja - porównanie x(t) oraz x_rec(t);
figure(1);
subplot(1,2,1);
plot(t,x,'r'); grid on; hold on;
plot(tn,xn,'bo');
xlim([-0.04,0.04]); xlabel('t [s]','fontsize',14); title('x(t) oraz x[n]','fontsize',14);
subplot(1,2,2);
plot(f,abs(X0),'r.-'); grid on; hold on;
plot(f,abs(X),'b-');
xlim([min(f),max(f)]); xlabel('f [Hz]','fontsize',14); title('Dt ^. abs(X(f)) oraz abs(X(F))','fontsize',14);
figure(2);
plot(t,x,'r.'); grid on; hold on;
plot(t,x_rec,'b');
xlim([-0.04,0.04]); xlabel('t [s]','fontsize',14); title('x(t) oraz x_r_e_c(t)','fontsize',14);
% KONIEC PRZYKŁADU 2.6.
Przykład 2.7
Przykład ten (plik przyklad_2_7.m) jest bardzo podobny do przykładu 2.6. Różnica polega na tym, że w tym przypadku pokazano jak zweryfikować oraz zilustrować pasmową wersję twierdzenia o próbkowaniu dla zadanego sygnału x(t). Sygnał pasmowy spełniający odpowiednie założenia otrzymano wykorzystując twierdzenie o modulacji - czyli mnożąc sygnał dolnopasmowy przez funkcję kosinus o częstotliwości fc. (w tym przypadku fc = 500 Hz).% Przykład 2.7:
% m-plik: przyklad_2_7.m
%
% Cyfrowe przetwarzanie sygnałów. Podstawy, multimedia, transmisja. PWN, Warszawa, 2014.
% Rozdział 2. Podstawy analizy częstotliwościowej i próbkowanie sygnałów.
% Autorzy rozdziału: Przemysław Korohoda, Krzysztof Duda.
%
% Więcej na stronie internetowej: http://teledsp.kt.agh.edu.pl
%
% Opracowanie przykładu: kwiecień 2013 r./PK.
clc; clear; close all;
dt=0.0001; tmax=0.5; t=-tmax:dt:tmax; % odcinek pseudo-ciągłej osi czasu;
t=t(:); % ułożenie elementów ciągu t w wektorze kolumnowym;
Dt=0.0025; fpr=1/Dt; % Dt to faktyczny odstęp próbkowania;
tn0=0:Dt:tmax; tn=[-tn0(end:-1:2),tn0]; % ciąg wartości t dla próbek, tak żeby w środku było t=0;
dF=0.0001; F=-2:dF:2; F=F(:); % odcinek pseudo-ciągłej osi F, F(:) - pionowo;
f=F*fpr; % odcinek pseudo-ciągłej osi częstotliwości f, w Hz;
f0=50; t0=0.01; M=2; % parametry sygnału;
fc=(fpr/2)*(2*M+1)/2; % fc - częstotliwość środkowa;
x=sinc(f0*(t-t0)).^2.*cos(2*pi*fc*t); % sygnał pseudo-ciągły: x(t);
xn=sinc(f0*(tn-t0)).^2.*cos(2*pi*fc*tn); % ciąg próbek: x[n];
NF=length(F); H=zeros(NF,1); % wstępne przygotowanie filtru odtwarzającego;
H(((f>-fpr/4-fc)&(f>fpr/4-fc))|((f>-fpr/4+fc)&(f>fpr/4+fc)))=1;
N=length(xn); X=zeros(NF,1); % przygotowanie do pętli liczącej DTFT;
for n=1:N,
n1=(n-(N+1)/2); % przeliczenie indeksów, by pasowały do wzoru DTFT;
X=X+xn(n)*exp(-j*2*pi*F*n1);
end;
X0=X.*H; % odtwarzamy widmo sygnału;
Nt=length(t); x_rec=zeros(Nt,1); % przygotowanie do pętli odtwarzającej sygnał;
for n=1:N,
n1=(n-(N+1)/2); % przeliczenie indeksów, by pasowały do wzoru;
t1=t-n1*Dt; % przesunięcie czasu dla każdego elementu sumy;
x_rec=x_rec+xn(n)*sinc(fpr/2*t1).*cos(pi*((2*M+1)/2)*fpr*t1);
end;
rec_err=max(abs(x_rec-x)), % weryfikacja - porównanie x(t) oraz x_rec(t);
figure(1);
subplot(1,2,1);
plot(t,x,'r'); grid on; hold on;
plot(tn,xn,'bo');
xlim([-0.04,0.04]); xlabel('t [s]','fontsize',14); title('x(t) oraz x[n]','fontsize',14);
subplot(1,2,2);
plot(f,abs(X0),'r.-'); grid on; hold on;
plot(f,abs(X),'b-');
xlim([min(f),max(f)]); xlabel('f [Hz]','fontsize',14); title('Dt ^. abs(X(f)) oraz abs(X(F))','fontsize',14);
figure(2);
plot(t,x,'r.'); grid on; hold on;
plot(t,x_rec,'b');
xlim([-0.04,0.04]); xlabel('t [s]','fontsize',14); title('x(t) oraz x_r_e_c(t)','fontsize',14);
% KONIEC PRZYKŁADU 2.7.
Przykład 2.8
Przykład 2.8 (plik przyklad_2_8.m) pokazuje jak: a) szybko wyznaczyć przybliżenie DTFT (poprzez uzupełnienie ciągu wartościami zerowymi i policzenie DFT), b) wykorzystać wynik DFT, by uzyskać opis osi poziomej w Hz (funkcja fftshift), c) zmieniać cechy fontów w opisie osi oraz wykresu, d) narzucać położenie wartości opisujących na osi wykresu, e) wykorzystać funkcję mod (modulo) w celu otrzymania próbek sygnału piłokształtnego, f) zapełniać wnętrze markerów oraz zmieniać grubość linii, g) opracować program dwuwariantowo – dla wartości N parzystej lub nieparzystej (warunek if).% Przykład 2.8:
% m-plik: przyklad_2_8.m
%
% Cyfrowe przetwarzanie sygnałów. Podstawy, multimedia, transmisja. PWN, Warszawa, 2014.
% Rozdział 2. Podstawy analizy częstotliwościowej i próbkowanie sygnałów.
% Autorzy rozdziału: Przemysław Korohoda, Krzysztof Duda.
%
% Więcej na stronie internetowej: http://teledsp.kt.agh.edu.pl
%
% Opracowanie przykładu: kwiecień 2013 r./PK.
clc; clear; close all;
N1=2^4-1;
x=mod(0:N1-1,5); % przykładowy ciąg próbek (warto zmienić);
X1=fft(x); % wyznaczamy DFT;
K=32; N2=K*N1; % zakładamy wydłużenie ciągu K+1 razy;
X2=fft(x,N2); % przybliżone wyznaczanie DTFT;
fpr=10; df1=fpr/N1; df2=fpr/N2; % potrzebne do opisu osi DFT oraz DTFT w Hz;
f1=0:df1:(N1-1)*df1; f2=0:df2:(N2-1)*df2; % wektory opisu osi DFT oraz DTFT w Hz;
sX1=fftshift(X1); % przestawienie elementów DFT;
if N1/2==round(N1/2), % dla N1 parzystych;
sf1=[f1(N1/2+1:end)-fpr,f1(1:N1/2)]; % opracowanie opisu osi poziomej w Hz;
else % dla N1 nieparzystych;
sf1=[f1((N1+1)/2+1:end)-fpr,f1(1:(N1+1)/2)]; % opracowanie opisu osi poziomej w Hz
end;
sX2=fftshift(X2); % przestawienie elementów przybliżenia DTFT;
if N2/2==round(N2/2),
sf2=[f2(N2/2+1:end)-fpr,f2(1:N2/2)]; % opracowanie opisu osi poziomej w Hz;
else
sf2=[f2((N2+1)/2+1:end)-fpr,f2(1:(N2+1)/2)]; % opracowanie opisu osi poziomej w Hz
end;
figure(1); % tę ilustrację warto rozciągnąć na cały ekran;
subplot(3,1,1);
stem(0:N1-1,x); xlim([0,N1-1]); grid on;
set(gca,'fontsize',12,'xtick',0:N1-1); xlabel('n','fontsize',18,'fontangle','italic');
subplot(3,1,2);
plot(f2,abs(X2),'b.-'); grid on; hold on;
plot(f1,abs(X1),'ro','linewidth',2,'markerface','r');
set(gca,'fontsize',12); xlabel('f [Hz]','fontsize',18,'fontangle','italic');
legend('DTFT','DFT');
subplot(3,1,3);
plot(sf2,abs(sX2),'b.-'); grid on; hold on;
plot(sf1,abs(sX1),'ro','linewidth',2,'markerface','r');
set(gca,'fontsize',12); xlabel('f [Hz]','fontsize',18,'fontangle','italic');
legend('DTFT','DFT');
% KONIEC PRZYKŁADU 2.8.
Przykład 2.9
Przykład (plik przyklad_2_9.m) pokazuje: a) jak wytworzyć i zbadać sygnał o zadanym przebiegu częstotliwości, b) czym rożni się informacja pozyskana z całościowego wyznaczania DFT od otrzymanej z DFT dla kolejnych bloków, c) co daje zastosowanie okna Hamminga, d) efekt aliasingu, e) jak wykorzystać wykresy trójwymiarowe do przedstawiania danych w macierzy, f) jak edytować opis wykresu.% Przykład 2.9:
% m-plik: przyklad_2_9.m
%
% Cyfrowe przetwarzanie sygnałów. Podstawy, multimedia, transmisja. PWN, Warszawa, 2014.
% Rozdział 2. Podstawy analizy częstotliwościowej i próbkowanie sygnałów.
% Autorzy rozdziału: Przemysław Korohoda, Krzysztof Duda.
%
% Więcej na stronie internetowej: http://teledsp.kt.agh.edu.pl
%
% Opracowanie przykładu: kwiecień 2013 r./PK.
clc; clear; close all;
N=1024; M=32; dM=4; fpr=100; % wartości parametrów;
Df=5; dt=1/fpr; f0=25; % wartości parametrów (cd.);
t=0:dt:(N-1)*dt; % oś czasu;
p=2*pi*(f0+0.5*Df*t).*t; % faza zależna od czasu;
x=cos(p); x=x(:); XN=fft(x); % sygnał i jego całościowe DFT;
w=hamming(M); k=0; % okno i inicjalizacja licznika;
for m0=1:dM:N-M+1; k=k+1;
x1=x(m0:m0+M-1); x2=x1.*w; % kolejne bloki i mnożenie przez okno;
X1(:,k)=fft(x1); X2(:,k)=fft(x2); % kolejne transformaty DFT;
end;
f=(0:M/2)*fpr/M; fN=(0:N/2)*fpr/N; n=(1:dM:N-M+1)+M/2; % przygotowanie dla wykresów;
fs_1=14; fs_2=12; % dwie wielkosci czcionek (font size);
figure(1); % tę ilustrację warto rozciągnąć, ale nie na cały ekran, tylko trochę;
% Uwaga - opisy osi wykresów 3D (mesh) mogą się pojawić trochę za nisko;
subplot(3,1,1);
plot(fN,abs(XN(1:N/2+1)),'ko'); grid on;
xlabel('f [Hz]','fontsize',fs_1); title('abs(X_N)','fontsize',fs_1);
set(gca,'fontsize',fs_2);
subplot(3,2,3);
mesh(t(n),f,abs(X1(1:M/2+1,:))); axis tight; view([-30,15]);
xlabel('t [s]','fontsize',fs_1); ylabel('f [Hz]','fontsize',fs_1);
title('abs(X_1)','fontsize',fs_1); set(gca,'fontsize',fs_2);
subplot(3,2,4);
mesh(t(n),f,abs(X2(1:M/2+1,:))); axis tight; view([-30,15]);
xlabel('t [s]','fontsize',fs_1); ylabel('f [Hz]','fontsize',fs_1);
title('abs(X_2)','fontsize',fs_1); set(gca,'fontsize',fs_2);
subplot(3,2,5);
mesh(t(n),f,abs(X1(1:M/2+1,:))); axis tight; view([0,90]);
xlabel('t [s]','fontsize',fs_1); ylabel('f [Hz]','fontsize',fs_1);
title('abs(X_1)','fontsize',fs_1); set(gca,'fontsize',fs_2);
subplot(3,2,6);
mesh(t(n),f,abs(X2(1:M/2+1,:))); axis tight; view([0,90]);
xlabel('t [s]','fontsize',fs_1); ylabel('f [Hz]','fontsize',fs_1);
title('abs(X_2)','fontsize',fs_1); set(gca,'fontsize',fs_2);
colormap(gray(256)); % wszystkie wykresy będą teraz "na szaro";
% KONIEC PRZYKŁADU 2.9.
Przykład 2.10
Implementacja algorytmu IpDFT BY1 (plik BY1_LC.m) z korekcją przecieku widmowego (Algorytm 2.1) w postaci funkcji. Argumentami wywołania funkcji są tłumiony sygnał sinusoidalny (2.74) oraz liczba iteracji korekcji przecieku widmowego NI. Funkcja zwraca estymowane wartości częstotliwości, tłumienia, amplitudy i fazy modelu (2.74). Dla NI = 0 używany jest tylko algorytm IpDFT BY1. Funkcję można używać również do analizy sygnałów nietłumionych, tj. dla d = 0, oraz do analizy sumy (tłumionych) sygnałów sinusoidalnych (w takim przypadku należy wybrać właściwe maksimum lokalne DFT). Przykładowe wykorzystanie tej funkcji podano w przykładzie 2.11.function [we, de, Ae, pe] = BY1_LC(x,NI);
% [we, de, Ae, pe] = BY1_LC(x,NI);
%
% Leakage correction for BY1 interpolated DFT
% x - input signal x=A*cos(w*n+p).*exp(-d*n)
% NI - number of iterations, for NI=0 BY1 without LC is computed
% Estimated values:
% we - frequency (rad)
% de - damping
% Ae - amplitude
% pe - phase (rad)
%
% Przykład 2.10.
% m-plik (funkcyjny): BY1_LC.m
%
% Cyfrowe przetwarzanie sygnałów. Podstawy, multimedia, transmisja. PWN, Warszawa, 2014.
% Rozdział 2. Podstawy analizy częstotliwościowej i próbkowanie sygnałów.
% Autorzy rozdziału: Przemysław Korohoda, Krzysztof Duda.
%
% Więcej na stronie internetowej: http://teledsp.kt.agh.edu.pl
%
% Opracowanie przykładu: kwiecień 2013 r./KD.
N = length(x);
Xw = fft(x); % computation of DFT
[Xabs, ind] = max(abs(Xw(1:round(N/2))));
k = [ind-1 ind ind+1];
dw = 2*pi/N; %DFT frequency step
wkm1= (k(1)-1)*dw; %frequency of DFT bin with index k-1
wk = (k(2)-1)*dw; %frequency of DFT bin with index k
wkp1= (k(3)-1)*dw; %frequency of DFT bin with index k+1
[we, de, Ae, pe, lam] = BY1_in_LC(wkm1,wk,wkp1,Xw(k),N);
%% Leakage correction for negative frequencies
wkk = [wkm1 wk wkp1];
for iter=1:NI
for m=1:3;
Xw_correction(m) = ...
(Ae/2)*( exp(-j*pe)*(1-conj(lam)^N)/(1-conj(lam)*exp(-j*wkk(m))));
end
Xw_correction = Xw(k) - Xw_correction; %
[we, de, Ae, pe, lam] = BY1_in_LC(wkm1, wk, wkp1, Xw_correction, N);
end
function [we, de, Ae, pe, lam] = BY1_in_LC(wkm1, wk, wkp1, Xw, N);
r = ( -exp(-j*wk)+exp(-j*wkm1) )/( -exp(-j*wkp1)+exp(-j*wk) );
R = ( Xw(1)-Xw(2) )/( Xw(2)-Xw(3) );
lam= exp(j*wk)*(r-R)/( r*exp(-j*2*pi/N)-R*exp(j*2*pi/N) );
we = imag(log(lam));
de = -real(log(lam));
if round(1e6*R)==-1e6 %%coherent sampling, d=0
Ae = 2*abs(Xw(2))/N;
pe = angle(Xw(2));
else
c = (1-lam^N)/(1-lam*exp(-j*wk));
c = 2*Xw(2)/c;
Ae = abs(c);
pe = angle(c);
end
% KONIEC PRZYKŁADU 2.10.
Przykład 2.11
Estymacja parametrów tłumionego sygnału sinusoidalnego (plik przyklad_2_11.m).% Przykład 2.11:
% m-plik: przyklad_2_11.m
%
% Cyfrowe przetwarzanie sygnałów. Podstawy, multimedia, transmisja. PWN, Warszawa, 2014.
% Rozdział 2. Podstawy analizy częstotliwościowej i próbkowanie sygnałów.
% Autorzy rozdziału: Przemysław Korohoda, Krzysztof Duda.
%
% Więcej na stronie internetowej: http://teledsp.kt.agh.edu.pl
%
% Opracowanie przykładu: kwiecień 2013 r./KD.
close all; clear all; clc;
fp= 100; % częstotliwość próbkowania w Hz;
N = 256; % liczba próbek;
A = 3; % amplituda sygnału;
f = 3; % częstotliwość sygnału analogowego w Hz;
p = 0.1; % faza w rad;
D = 0.1; % tłumienie;
dt= 1/(fp);
t = (0:N-1)*dt;
%% sygnał testowy:
x = A*cos(2*pi*f*t+p).*exp(-D*t);
%% BY1 bez korekcji przecieku widmowego:
NI=0; [Ome, de, Ae, pe] = BY1_LC(x,NI);
fe = Ome*fp/(2*pi); % przeliczenie modelu;
De = de*fp; % przeliczenie modelu;
blad(1,:) = [fe-f, De-D, Ae-A, pe-p];
%% BY1 korekcja przecieku widmowego 5 iteracji:
NI=5; [Ome, de, Ae, pe] = BY1_LC(x,NI);
fe = Ome*fp/(2*pi); % przeliczenie modelu;
De = de*fp; % przeliczenie modelu;
blad(2,:) = [fe-f, De-D, Ae-A, pe-p];
%%
blad_f_D_A_p = num2str(blad),
%% Rysunki
%% DFT
Xd = fft(x);
[val ind] = max(abs(Xd(1:round(N/2))));
dw = (2*pi/N);
wd = (0:N-1)*dw;
fd = fp*wd/(2*pi);
%% DTFT
Nfft = 2^16;
Xc = fft(x,Nfft); % zagęszczone próbkowanie widma sygnału dyskretnego przez dołaczenie wektora zer do długości sygnału Nfft przed liczeniem DFT
[val ind] = max(abs(Xc(1:round(Nfft/2))));
dwc= (2*pi/Nfft);
wdc= (0:Nfft-1)*dwc;
fdc= fp*wdc/(2*pi);
FontSize=16;
FontName = 'Times New Roman';
figure, hold on
plot(t,x,'.-k','LineWidth',1,'MarkerSize',12),
title('x(t)= 3cos(2\pi3t+0.1)exp(-0.1t)','FontSize',FontSize,'FontName',FontName)
xlabel('t_n [s]','FontSize',FontSize,'FontName',FontName)
ylabel('x(t_n)','FontSize',FontSize,'FontName',FontName)
axis tight,box on
v=axis; axis([v(1) v(2) -3.01 3.01])
set(gca,'FontSize',FontSize,'FontName',FontName)
figure, hold on
H1=stem(fd, abs(Xd),'ob','filled'); set(H1,'LineWidth',1);
plot(fdc,abs(Xc),'-r','LineWidth',1),
H2=stem(fe, 400, '--k','filled'); set(H2,'LineWidth',2);
plot(fdc,abs(Xc),'-r','LineWidth',1),
legend('|DFT|','|DTFT|')
xlabel('f [Hz]','FontSize',FontSize,'FontName',FontName)
axis tight,box on
set(gca,'FontSize',FontSize,'FontName',FontName)
v=axis; axis([2 6 v(3) 350])
% KONIEC PRZYKŁADU 2.11.
Przykład 2.12
W przykładzie (plik przyklad_2_12.m) zademonstrowano weryfikację wyprowadzonych zależności na opóźnienie grupowe i fazowe. Zmieniając wybrane parametry, można zbadać efekty wynikające z założeń upraszczających. Odpowiednią ilustrację graficzną należy zrealizować samodzielnie. Warto obejrzeć charakterystyki kanału, przebiegi czasowe sygnałów oraz różnicy między wynikiem filtrowania i sygnału otrzymanego po zastosowaniu wyliczonych opóźnień. Należy zwrócić uwagę, że przykład pokazuje obliczenia wykonane w dziedzinie sygnałów dyskretnych.% Przykład 2.12:
% m-plik: przyklad_2_12.m
%
% Cyfrowe przetwarzanie sygnałów. Podstawy, multimedia, transmisja. PWN, Warszawa, 2014.
% Rozdział 2. Podstawy analizy częstotliwościowej i próbkowanie sygnałów.
% Autorzy rozdziału: Przemysław Korohoda, Krzysztof Duda.
%
% Więcej na stronie internetowej: http://teledsp.kt.agh.edu.pl
%
% Opracowanie przykładu: kwiecień 2013 r./PK.
close all; clear all; clc;
N=2^13; fpr=1000; % długość ciągu i częstotliwość próbkowania w [Hz];
dt=1/fpr; t=0:dt:(N-1)*dt; % oś czasu w [s];
[b1,a1]=butter(15,0.6); % modelujemy kanał jako filtr Butterwortha (patrz rozdział 4);
[H,W]=freqz(b1,a1,5000); % wyznaczamy charakterystyki kanału transmisyjnego;
F=W/(2*pi); f=fpr*F; % przeliczamy częstotliwości;
%% definiujemy sygnał modulujący (xm) oraz nośną (xc):
fm=2; fc=300; % przyjmujemy częstotliwości sygnału modulującego i nośnej w [Hz];
xm=sinc(fm*(t-(N/2)*dt)); % sygnał modulujący;
xc=cos(2*pi*fc*t); % sygnał nośnej;
x=xm.*xc; % sygnał zmodulowany;
%% "przesyłamy" sygnał przez modelowany cyfrowo kanał transmisyjny:
y=filter(b1,a1,x);
P=unwrap(angle(H)); df=f(2)-f(1); % wyliczamy przesunięcie fazy i przyrost częstotliwości;
k1=find((f>fc-df/2)&(f<=fc+df/2)); % wyznaczamy indeks dla f=fc;
A=abs(H(k1)); % zmiana amplitudy sygnału dla f=fc;
tp=-P(k1)/(2*pi*fc), % opóźnienie fazowe w [s] według wzoru;
tg=-(1/(2*pi))*(P(k1+1)-P(k1-1))/(2*df), % opóźnienie grupowe w [s] według wzoru;
tg_vs_dt=tg/dt, % proporcja tg do kroku próbkowania na osi czasu;
tp_vs_dt=tp/dt, % proporcja tp do kroku próbkowania na osi czasu;
%% na podstawie wyliczonych wartości opóźnień konstruujemy sygnał porównawczy:
xm_test=sinc(fm*(t-(N/2)*dt-tg)); % opóźniony sygnał modulujący;
xc_test=cos(2*pi*fc*(t-tp)); % opóźniony sygnał nośnej;
y_test=A*xc_test.*xm_test; % sygnał zmodulowany po transmisji;
x_err=y_test-y; % sygnał różnicy (błędu);
max_abs_err=max(abs(x_err)), % maksymalny błąd bezwzględny;
figure(1); % fragment charakterystyki amplitudowej filtra modelującego kanał transmisyjny;
k=k1-20:k1+20;
plot(f(k),abs(H(k)),'b.-'); grid on; hold on; axis tight;
plot([f(k(1)),f(k(end))],[abs(H(k(1))),abs(H(k(end)))],'k');
xlabel('f [Hz]'); title('Fragment ch-ki amplitudowej kanalu w rejonie f nosnej');
figure(2); % graficzne porównanie sygnałów;
plot(t,y,'b.-'); grid on; hold on;
plot(t,y_test,'r.');
legend('y(t)','y_t_e_s_t(t)'); xlabel('t [s]'); title('Graficzne porownanie sygnalow (prosze uzyc "zoom")');
axis tight;
% sygnały są bardzo podobne, ale nie identyczne - dlaczego?
% KONIEC PRZYKŁADU 2.12.