Laboratorul 04.

Semnale în domeniul frecvență

Materiale ajutătoare:

    1. Secțiunile 4.5 (Exercițiul 1), 4.7 (Exercițiul 2), 4.6 (Exercițiul 3)

Exercițiul 1 - aproximare de semnale [4p]

O să folosim din nou semnalul dreptunghiular, de amplitudine 'A' în intervalul [0, T/2] și '-A' în intervalul [T/2, T] care are coeficienții Fourier dați de formula:

\begin{equation} c_{k} = \left\lbrace \begin{array}{} \frac{2}{j \pi k}A \qquad k \quad impar \\0 \qquad \quad k \quad par \end{array} \right. \end{equation}

În acest exercițiu vom încerca să calculăm rădăcina pătratică medie (eng. root mean square) a semnalului de eroare $\epsilon_{K}$ dat de: \begin{equation} \epsilon_{K}(t) = s(t) - s_{K}(t), \end{equation}

unde $s_{K}(t)$ este semnalul aproximat cu K termeni. Folosind transformata Fourier avem: \begin{equation} s_K(t) = \sum_{k=-K}^K c_k e^\frac{j 2 \pi k t}{T} \end{equation}

iar RMS-ul lui $\epsilon_K$ este dat de: \begin{equation} \text{rms}(\epsilon_K) = \sqrt{2\cdot \sum_{k=K+1}^\infty |c_k|^2} \end{equation}

Task-ul vostru este să determinați valoarea lui K astfel încât $\text{rms}(\epsilon_K)$ este aproape 0 și după aceea să vedeți că într-adevăr semnalul reconstruit aproximează bine semnalul original.

Pentru asta ar trebui să urmăriți următorii pași:

  1. Creați semnalul original. Folosiți, de exemplu T=100 și A=1 pentru a genera semnalul cu valoarea 1 în primele 50 de eșantioane și -1 în ultimele 50. Reprezentaţi grafic semnalul ca funcţie de timp (unde timpul începe de la 1 până la 100) [1p]. Puteți să vă folosiți de acest cod:
    A = 1;
    T = 100;
    x = 1:T;
    s = -A*ones(1, T);
    s(1:(T/2)) = A;
  2. Calculați coeficienții Fourier $c_{k}$ pentru $k=[-10:10]$. Plotați amplitudinile $|c_k|^2$ ca funcție de $k$. [1p]
  3. Calculați coeficienții Fourier $c_{k}$ pentru $k=[0:1500]$. Calculați $\text{rms}(\epsilon_{K})$ ca funcție de $K$ pentru $K \in \{1, \ldots, 500\}$. Plotați valoarea rms pentru $K \in \{1, \ldots, 500\}$.[1p]
  4. Determinați cel mai mic $K$ astfel încât $\text{rms}(\epsilon_{K}) < 0.025$ și reconstruiți semnalul original folosind acest număr de coeficienți. Trebuie să folosiți atât coeficienții pozitivi cât și negativi (de ex. de la -K la K) pentru a reconstrui semnalul. Reprezentați grafic semnalul reconstruit și comparați-l cu semnalul inițial. [1p]

Semnalele reale au următoarea proprietate: coeficienții Fourier negativi sunt conjugații complexi ai celor pozitivi $c_{-k} = c_{k}^*$. De asemenea, semnalele impare s(-t) = -s(t), au coeficienți complet imaginari, obținând $c_{-k} = -c_{k}$

Pentru a calcula rms al erorii trebuie sa calculam suma pentru toti coeficienții $c_k$ cu $|k|> k_0 $, adica o infinitate de termeni. Putem încerca doar să aproximăm această sumă, sau ne putem folosi de unele proprietăți ale serier Fourier pentru a o calcula exact. Mai precis, vom folosi Teorema lui Parseval prin care putem calcula puterea unui semnal in doua feluri, in domeniul timp, integrand semnalul la pătrat peste o perioadă sau în frecvență calculand suma pătratelor modulului ale fiecărui coeficient

\begin{equation} \frac{1}{T}\int_0^T{s^2(t)dt} = \sum_{-\infty}^{\infty}{|c_k|^2} \end{equation}

Folosind cele descrise mai sus puteți calcula exact rms-ul erorii:

  1. calculați puterea totala 'de mână' pentru semnalul dreptunghiular.
  2. prin scaderea pătratelor termenilor de la $c_k, k \in \{-k...k\}$

O posibilă soluție:

ex1.m
A = 1;
T = 100;
x = 1:T;
s = -A*ones(1, T);
s(1:(T/2)) = A;
 
h = figure;
plot(x,s);
ylabel('Amplitude');
xlabel('Time');
ylim([-A-0.5, A+0.5]);
title('Original signal s(t)');
print(h, '-dpng', 'original_signal.png');
 
% Compute some of the Fourier coefficients
cmax = 21;
kv = -cmax:cmax;
N = length(kv);
coef = zeros(N,1);
for i=1:N
    k = kv(i);
    if mod(k,2) ~= 0
        coef(i) = (2*A) / (j*pi*k);
    end
end
 
% Plot power of coefficients
h = figure;
csquare = coef.*coef;
stem(kv,abs(csquare));
xlabel('Frequency component (k)');
ylabel('Magnitude of component');
title('Fourier coefficients squared');
print(h, '-dpng', 'squared_coefficients.png');
 
% Plot rms of error signal
cmax = 1000;
kv = 0:cmax;
N = length(kv);
coef = zeros(N,1);
for i=1:N
    k = kv(i);
    if mod(k,2) ~= 0
        coef(i) = (2*A) / (j*pi*k);
    end
end
 
kmax = 500;
rmse = zeros(kmax, 1);
for i=1:kmax
    ss = coef(i:end)'*coef(i:end);
    rmse(i) = sqrt(2*ss);
end
 
h=figure;
plot(0:(kmax-1), rmse); % Look with plot, semilogy and loglog
title('RMS of e_k');
ylabel('RMS(e_k)');
xlabel('K');
grid on;
print(h, '-dpng', 'rms_ek.png');
 
% Reconstruct signal with only some of the components
% Note: need to use both positive and negative k's
% The fourier coeficients  c(-k) = conj(c(k)), for real signals
% As the original signal is odd, the fourier coeficient are completly imaginary
% so c(-k) = -c(k)
 
fk = find(rmse < 0.03);
N = fk(1);
sr = zeros(1, T);
for t=1:T
    for i=1:N
        sr(t) = sr(t) + coef(i)*exp(j*2*pi*kv(i)*t/T) ...
                      - coef(i)*exp(-j*2*pi*kv(i)*t/T);
    end
end
 
h=figure;
plot(1:T, sr);
ylim([-A-1, A+1]);
title('Reconstructed signal s(t)');
ylabel('Amplitude');
xlabel('Time');
print(h, '-dpng', 'reconstructed_signal.png');

Exercițiul 2 - filtrare [4p]

În acest exercițiu trebuie să calculăm coeficienții Fourier ai output-ului unui filtru trece-jos, dat fiind un semnal de intrare de tip puls.

Pentru un semnal dat putem găsi coeficienții Fourier ($c_k$). Atunci, dacă știm funcția de transfer a filtrului trece-jos (sau alt tip de sistem liniar), am arătat la curs că putem găsi coeficienții Fourier ($c_k^y$) ai semnalului rezultat ca: \begin{equation} c_k^y = H(\frac{k}{T}) \cdot c_k \end{equation}

Astfel putem reconstrui semnalul de ieșire folosind coeficienții Fourier $c_k^y$.

Ni se dă ca input un semnal de tip puls cu amplitudine $A=1$ și pulsul de durată $\Delta=\frac{T}{5}$. Știm că coeficienții Fourier ai semnalului sunt dați de: \begin{equation} c_k = A \cdot e^{-j\frac{\pi k \Delta}{T}} \cdot \frac{\text{sin}(\frac{\pi k \Delta}{T})}{\pi k} = A \cdot e^{-j\frac{\pi k \Delta}{T}} \cdot \frac{\Delta}{T} \cdot \text{sinc}(\frac{\pi k \Delta}{T}). \end{equation}

Atenție: Matlab folosește funcția sinc normalizată sinc(x) = $ \frac{\text{sin}(\pi x)}{\pi x}$

Puteți vedea circuitul filtrului trece-jos în următoarea imagine:

Funcția de transfer a circuitului (pe care am determinat-o la curs și la laboratoarele anterioare) este următoarea: \begin{equation} H(f=\frac{k}{T}) = \frac{1}{1+j 2 \pi R C \frac{k}{T}} \end{equation} unde R și C sunt rezistența și respectiv capacitatea.

Task-ul vostru este să determinați coeficienții output-ului și să reconstruiți semnalul de ieșire pentru diferite frecvențe de tăiere.

Pentru aceasta urmăriți următorii pași:

  1. Generați semnalul puls și plotați-l. Puteți folosi T=100 de eșantioane, dintre care doar $\Delta$ nu sunt egale cu 0. [1p]
  2. Calculați primii N=30 coeficienți Fourier pozitivi $c_k$ ai semnalului și plotați-i. [1p]
  3. Calculați coeficienții Fourier asociați semnalului de output, $c_k^y$, folosind formula de mai sus și plotați-i. Pentru asta va trebui să alegeți o frecvență de cut-off $f_c$ care va determina valorile R și C ( $RC = \frac{1}{2\pi f_c}$ ). [1p]

Puteți încerca următoarele valori pentru $f_c$:

  1. $f_c = 0.1/T$ (frecvența de cut-off e mult mai mică decât frecvența fundamentală a semnalului ⇒ filtrare puternică)
  2. $f_c = 1/T$ (frecvența de cut-off = frecvența fundamentală ⇒ puterea este înjumătățită)
  3. $f_c = 10/T$ (frecvența de cut-off mult mai mare decât frecvența fundamentală ⇒ filtrare slabă)
  4. Reconstruiți semnalul de output cu ajutorul seriei Fourier (folosind formula de la exercițiul 1) [1p]

Ideea este să folosiți seriile Fourier pentru a determina ieșirea filtrului linear.

ex3.m
A = 1;
D = 20;
T = 100;
x = 1:T;
s = zeros(1,T);
s(1:D) = A;
 
h = figure;
plot(x,s);
ylabel('Amplitude');
xlabel('Time');
ylim([-A-0.5, A+0.5]);
title('Original signal s(t)');
print(h, '-dpng', 'pulse_signal.png');
 
% Compute some of the Fourier coefficients
cmax = 20;
kv = 0:cmax;
N = length(kv);
coef = zeros(N,1);
for i=1:N
    k = kv(i);
    coef(i) = (D/T)*A*exp(-j*pi*k*D/T)*sinc(k*D/T);
end
 
% Plot magnitude of pulse's coefficients
h = figure;
stem(kv, abs(coef));
xlabel('Frequency component (k)');
ylabel('Magnitude of component');
title('Input Fourier coefficients');
hold on;
plot(kv, abs(coef), 'r--');
print(h, '-dpng', 'coef_pulse.png');
 
% Compute coefficients of output after lp filter
fcoef = zeros(N,1);
%fc = 0.1/T;
%fc = 1/T;
fc = 10/T;
RC = 1/(2*pi*fc);
for i=1:N
    k = kv(i);
    fcoef(i) = coef(i) / (1+j*2*pi*RC*k/T);
end
 
% Plot magnitude of output's coefficients
h = figure;
stem(kv, abs(fcoef));
xlabel('Frequency component (k)');
ylabel('Magnitude of component');
title('Output Fourier coefficients');
hold on;
plot(kv, abs(fcoef), 'r--');
print(h, '-dpng', 'coef_out.png');
 
%Construct output signal with the computed coefficients
sr = zeros(1, T);
for t=1:T
    for i=1:N
        if kv(i) ~= 0
            sr(t) = sr(t) + fcoef(i)*exp(j*2*pi*kv(i)*t/T) ...
                    + conj(fcoef(i))*exp(-j*2*pi*kv(i)*t/T);
        else
            sr(t) = sr(t) + fcoef(i)*exp(j*2*pi*kv(i)*t/T);
        end
    end
end
 
h=figure;
plot(1:T, sr);
ylim([-A, A]);
title('Output filtered signal s(t)');
ylabel('Amplitude');
xlabel('Time');
print(h, '-dpng', 'output_signal.png');

If time allows and some of the students are interested, it would be nice to make a SPICE simulation (e.g. on ngspice.com) with a pulse signal and an RC low-pass filter to check that the output of the simulation is similar to what we did in MATLAB.

Here is an example of what to write in ngspice to get the “shark teeth” (for fc=ft):

*RC Filter Transient example
	
V1 in 0 PULSE(0 1 1ns 1ns 1ns 40us 200us)
R1 in out 3.2k
C1 out 0 10n
	
.TRAN 10n 1m
.end

Change the capacitor value to modify the effect of the filter (smaller values mean a weaker filtering, while larger values mean a stronger filtering).

Exercițiul 3 - comunicație digitală [2p]

Am văzut la curs că pentru a transmite 2 biți simultan putem folosi două frecvențe diferite (f1, f2) pentru a coda o valoare de 2 biți:

  1. '00': folosim un semnal egal cu 0 (nicio frecvență)
  2. '01': folosim o sinusoidă ce conține doar prima frecvență ($\sin(2\pi f_1 t)$)
  3. '10': folosim o sinusoidă ce conține doar a doua frecvență ($\sin(2\pi f_2 t)$)
  4. '11': folosim ambele frecvențe $f_1$ și $f_2$

Task-ul vostru e să creați o secvență random de 10 valori între 0 și 3 (pentru a folosi toate valorile de mai sus) și apoi să o codificați folosind 2 sinusoide așa cum este descris mai sus. Pentru asta ar trebui să:

  1. selectați frecvențele $f_1$ și $f_2$ astfel încât ele folosesc aceeași frecvență fundamentală (de ex.: $f_1 = 1 \cdot f_t$, $f_2 = 2 \cdot f_t$);
  2. plotați semnalul rezultat folosind perioada ($1/f_t$) pentru fiecare valoare transmisă;
  3. verificați că semnalul rezultat codează secvența voastră random;

Notă: pentru a genera o secvență random de valori întregi inspectați funcția 'randi' din MATLAB.

Posibila solutie:

%% Frquency modulation using two frequencies
close all;
k1 = 1;
k2 = 2;
T = 1;
f1 = k1/T;
f2 = k2/T;
t = 0:T/100:T;
s1 = sin(2*pi*f1*t);
s2 = sin(2*pi*f2*t);
values = randi([0 3], 1, 10);
nv = length(values);
s = [];
for i=1:nv
    si = zeros(1,length(s1));
    if mod(values(i), 2) == 1
        si = si + s1;
    end
    if values(i) > 1
        si = si + s2;
    end
    s = [s, si];
end
tt = length(s);
x = ((1:tt) - 1) * nv / tt;
h1 = figure;
plot(x, s);
grid on;
xlabel('Time (s)');
ylabel('Amplitude');
title('Frequency modulated signal');
print(h1, '-dpng', 'freqmod.png');
fprintf('The transmitted values are:\n');
values

ps/laboratoare/04.txt · Last modified: 2018/11/05 13:34 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