Fondamenti di Segnali e Trasmissione
Terzo laboratorio
Trasformata di Fourier di segnali continui in MATLAB
La funzione MATLAB che calcola la sequenza che approssima la trasformata di Fourier di un segnale è:
X=fft(x)*dt;
Il termine moltiplicativo dt permette di avere ampiezze congruenti con la definizione di trasformata di
Fourier di segnali continui. La lunghezza del vettore X, trasformata di Fourier, è uguale alla lunghezza del
segnale nei tempi (sequenza x). L’operazione di trasformazione di Fourier inversa è invece:
x=ifft(X)/dt;
I segnali discreti (sequenze nel tempo con passo di discretizzazione dt) che possono essere rappresentate in
MATLAB, hanno trasformate periodiche in frequenza, di periodo 1/dt. La sequenza X, output della
funzione fft() rappresenta l’osservazione su un solo periodo della trasformata di Fourier. Le frequenze su
cui è definita la sequenza X vanno quindi da –1/(2*dt) a +1/(2*dt). Il primo estremo è incluso solo se la
lunghezza del vettore x è un numero pari, mentre entrambi gli estremi sono esclusi per lunghezze dispari.
Le convenzioni adottate da MATLAB sono che il primo campione di x, x(1), è quello che corrisponde a t=0,
ed il primo campione di X, X(1) è quello che corrisponde a f=0. Esiste però la funzione fftshift() che
permette di riordinare le frequenze.
Quindi per poter calcolare la trasformata di Fourier del segnale x(t) (approssimato dalla sequenza x):
» N=length(x);
» X=fftshift(fft(x))*dt;
» df=1/(N*dt);
» f=[-N/2+[0:N-1]]*df; % se N è pari
» f=[-(N-1)/2+[0:N-1]]*df; % se N è dispari
Per visualizzare il risultato è importante ricordare che la trasformata X sarà sempre un vettore complesso.
Quindi si potranno mostrare modulo e fase:
» subplot(2,1,1), plot(f,abs(X));
» subplot(2,1,2), plot(f,angle(X));
oppure parte reale e parte immaginaria:
» subplot(2,1,1), plot(f,real(X));
» subplot(2,1,2), plot(f,imag(X));
Esercizio 0: La TDF del rettangolo ritardato
Si consideri il segnale rettangolare tempo continuo:
s(t)=A*rect( (t-T/2)/T )
dove A=3 e T=0.1 s. Si generi tale segnale tramite Matlab, ricorrendo alla funzione rect() definita nelle
lezioni precedenti. Di tale segnale, si calcoli quindi la Trasformata di Fourier e se ne visualizzi il
modulo e la fase. Si confronti, infine, quanto ottenuto con il risultato teorico.
Come asse temporale si definisca il vettore t=[0:0.001:4], espresso in secondi.
Per la risoluzione dell’esercizio, viene proposto il seguente codice Matlab:
dt=.001; % passo di discretizzazione asse temporale
T=0.1; % durata del rettangolo
A =3; % ampiezza del rettangolo
t =[0:dt:4];
R =A*rect(t,0,T-dt); % il rettangolo nei tempi
Rf=fftshift(fft(R))*dt;
N =length(t);
df=1/(N*dt);
f=[-(N-1)/2+[0:N-1]]*df; % asse delle frequenze
RfT=(A*T*sin(pi*f*T)./(pi*f*T)).*exp(-j*2*pi*f*((T-dt)/2)); % TDF teorica
RfT((length(f)+1)/2)=A*T; % MATLAB non risolve limiti notevoli!!
figure,
subplot(2,1,1),
plot(f,abs(RfT),'-b'),
hold on,
plot(f,abs(Rf),'.r'),
title('Modulo')
subplot(2,1,2),
plot(f,angle(RfT),'-b'),
hold on,
plot(f,angle(Rf),'.r'),
title('Fase')
xlabel('Frequenza [Hz]');
Modulo
0.4
0.3
0.2
0.1
0
-500
-400
-300
-200
-100
0
100
200
300
400
500
-100
0
100
Frequenza [Hz]
200
300
400
500
Fase
4
2
0
-2
-4
-500
-400
-300
-200
Esercizio 1: L’onda quadra
Note teoriche. Si consideri un’onda quadra di periodo 2T con duty cycle = 50% e media 1/2. Tale
segnale può essere generato replicando a passo 2T il segnale rect(t/T). Dalla teoria, è noto che la
trasformata di Fourier del segnale periodico così generato si ottiene campionando a passo 1/2T e
opportunamente scalando la trasformata di Fourier del segnale originale.
Quanto detto è espresso meglio nel grafico seguente:
T
-2T
-T
T
2T
media
−
−
5
2T
3
2T
3
2T
−
1
2T
1
2T
5
2T
Come si può notare, la sua trasformata di Fourier è composta da un numero infinito di armoniche
dispari, ossia di posizione in frequenza pari a 1//2T), 3/(2T), 5/(2T), …
Si produca tramite Matlab un codice per la generazione di un’onda quadra a media nulla, duty cycle
50% e periodo 2 s. Si studi in primo luogo la risposta in frequenza di tale segnale. Cosa succederebbe
se si riducesse il numero di armoniche della risposta in frequenza?
Per produrre l’onda quadra desiderata, si propone la soluzione seguente. Si scrivano le funzioni Matlab:
function y=even(x)
L=length(x(:)); % lunghezza del vettore
y=zeros(size(x));
for k=1:L,
if (floor(x(k)/2)==x(k)/2),
y(k)=1;
end
end
function y=ondaquadra(x,T), % T è il semiperiodo
y=-ones(size(x)); % ones(r,c) crea una matrice inizializzata a 1
set=find(even(floor(x/T))==1);
y(set)=1;
Si generi quindi l’onda quadra con i comandi:
dt=0.001;
t=[0:dt:8];
y=ondaquadra(t,1);
Onda quadra
1
0.8
0.6
0.4
0.2
0
-0.2
-0.4
-0.6
-0.8
-1
0
1
2
3
4
tempo [s]
5
Si visualizzi quindi la TDF di y:
Y=fftshift(fft(y))*dt;
N=length(t);
df=1/(N*dt);
asse_frequenze=[-(N-1)/2+[0:N-1]]*df; % N dispari
figure
plot(asse_frequenze,abs(Y));
title('Modulo della TDF')
xlabel('Frequenze [Hz]')
6
7
8
Modulo della TDF
5
4.5
4
3.5
3
2.5
2
1.5
1
0.5
0
0
1
2
3
4
5
6
Frequenze [Hz]
7
8
9
10
Per ridurre il numero di armoniche, per esempio alle prime 5, si utilizzi il seguente codice:
filter=[ones(1,41),zeros(1,length(Y)-81),ones(1,40)]; % prime 5 armoniche
Y_filtered=fft(y)*dt.*filter;
y_filtered=ifft(Y_filtered)/dt;
figure(4)
plot(t,y,'r') % onda quadra originale
hold on
plot(t,y_filtered) % onda quadra filtrata
Si utilizzi poi il codice appena scritto per filtrare diversi numeri di armoniche.
Onda quadra con prime 5 armoniche
1.5
1
0.5
0
-0.5
-1
-1.5
0
1
2
3
4
tempo [s]
5
6
7
8
7
8
Onda quadra con prime 10 armoniche
1.5
1
0.5
0
-0.5
-1
-1.5
0
1
2
3
4
tempo [s]
5
6
Onda quadra con prime 30 armoniche
1.5
1
0.5
0
-0.5
-1
-1.5
0
1
2
3
4
tempo [s]
5
6
7
8
Esercizio 2: Chiamata telefonica: Dual-Tone Multi Frequency
Le cifre del telefono a toni sono codificate con una coppia di sinusoidi. Premendo un tasto viene
generato il segnale:
y(t)=cos(2πfrt)+cos(2πfct)+w(t)
dove si assume la presenza di un disturbo casuale gaussiano a media nulla con varianza σ2. Inoltre il
tasto premuto è identificato dal valore delle due frequenze secondo lo schema seguente:
Si assuma che il segnale a tempo continuo venga campionato con la frequenza f0 = 8kHz, e che si
osservi il segnale per un periodo di osservazione Toss= 128 ms (corrispondente al tempo di pressione
del tasto).
Per identificare il tasto premuto a partire dal segnale y(t), se ne può valutare la risposta nel dominio
della frequenza e selezionare i massimi (due in questo caso) che corrispondono alle due sinusoidi
Come applicazione, si vuole simulare in Matlab la pressione di un dato tasto. Si consideri ad esempio il
tasto 4. Dai parametri precedenti la sequenza ha N=128ms*8kHz=1024 campioni. Il codice seguente
simula la pressione del tasto 4:
f0=8000; sigma2=2;
dt=1/f0;
fr=[697 770 852 941]; fc=[1209 1336 1477 1633];
frn=fr/f0; fcn=fc/f0;
r=2; c=1; % Tasto 4
N=1024; n=[0:N-1];
w=sqrt(sigma2)*randn(1,N); % Simulazione del rumore
y1=sin(2*pi*frn(r)*n) % simulazione della prima sinusoide
y2=sin(2*pi*fcn(c)*n) % simulazione della seconda sinusoide
y=y1+y2+w;
sound(y,f0); % Generazione del suono (si provi anche con y1 e y2)
Si determini visualizzi ora il contenuto in frequenza del segnale y(t), e se ne valuti il contenuto nel caso
senza rumore e con rumore (gaussiano a media nulla e varianza due)
dt=1/f0;
df=1/(N*dt);
f=[-(N)/2+[0:N-1]]*df; % N è pari
Y=fftshift(fft(y))*dt;
Modulo della risposta in frequenza
0.06
0.05
0.04
0.03
0.02
0.01
0
-4000
-3000
-2000
-1000
0
1000
frequenze [Hz]
2000
3000
4000
Si carichi ora il file y_es3 tramite il seguente comando:
load y_es3
A seguito di questa operazione, nel Workspace di Matlab viene caricato il vettore y. Esso corrisponde
alla pressione di un tasto incognito (come prima, f0=8000Hz, N=1024 campioni). Si applichi quanto
detto fino ad ora per determinare il tasto premuto.
Esercizio 3: Effetto Doppler
Lo scopo di questo esercizio è di valutare con un esempio pratico l’effetto Doppler.
Caricate il file siren.mat con il comando
load siren
A seguito di tale comando, nel workspace vengono caricate due variabili: sir1 e sir2. Tali sequenze
possono essere ascoltate con i comandi
fc=44100;
sound(sir1,fc)
sound(sir2,fc)
dove fc=44100Hz è la frequenza di campionamento di segnali audio di qualità CD. Come potete sentire,
si tratta di una sirena che si sta avvicinando all’ascoltatore (sir1) per poi allontanarsi (sir2). Per effetto
Doppler, mentre la sirena si avvicina la frequenza del segnale percepito dall’ascoltatore cresce mentre
decresce quando la sirena si allontana. Si verifichi quanto detto eseguendo un’analisi spettrale del
segnale tramite TDF. Il risultato è riportato nella figura seguente (legenda: in rosso il modulo della
TDF di sir2 e in blu il modulo della TDF di sir1).
Spunto di soluzione:
dt=1/fc;
N=length(sir1);
df=1/(N*dt);
f=[-(N/2)+[0:(N-1)]]*df;
S1=fftshift(fft(sir1))*dt;
S2=fftshift(fft(sir2))*dt;
figure(5)
plot(f,abs(S1))
hold on
plot(f,abs(S2),'r')
-3
2
Effetto Doppler
x 10
1.8
1.6
1.4
1.2
1
0.8
0.6
0.4
0.2
0
-6000
-4000
-2000
0
frequenza [Hz]
2000
4000
6000
Esercizio 4: Handel's Hallelujah Chorus
Si carichi e si ascolti il brano musicale Handel's Hallelujah Chorus (presente nelle librerie di Matlab)
tramite il comando:
fc=8192;
load handel
sound(y,fc)
Si visualizzi il modulo della trasformata di Fourier. Si filtri poi la sequenza musicale con un filtro a
media mobile del tipo
filt=ones(1,20)/20.
Si ascolti il segnale risultato dal filtraggio e se ne visualizzi il comportamento in frequenza,
confrontandolo con il precedente.
Modulo delle TDF
0.09
0.08
0.07
0.06
0.05
0.04
0.03
0.02
0.01
0
-5000 -4000 -3000 -2000 -1000
0
1000
frequenza [Hz]
2000
3000
4000
5000
© Copyright 2024 ExpyDoc
