Problema con Octave

[fairyvilje]

Sto lavorando ad una tesina per segnali e sistemi. Mi scuso per la scarsa qualità del codice .

Codice: Seleziona tutto

clear all
close all
clc

mf_c=32*10^3;

my=wavread("./Dati/segnaleDistorto_1.wav");
mt=1:length(my)

% Draw the signal in time
figure
plot(mt*1/mf_c,my)

grid
legend('Raw Wav')
axis([0 length(mt)*1/mf_c -1 1])


% Draw the signal in frequency
my_fft=fft(my);
my_current=abs(my_fft)*1/mf_c;
mf=1:length(my_current)

figure
plot(mf/length(my_current)*32000,my_current)

grid
legend('FFT')
axis([0 32000 0 max(my_current)])


% Draw the signal*cos(2*pi*f_o) in time
figure

for c=1:length(my)
  my2(c)=my(c)*cos(2*pi*1/mf_c*c*7500);
  c++
endfor

plot(mt*1/mf_c,my2);

grid
legend('Raw Wav')
axis([0 length(mt)*1/mf_c -1 1])


% Draw the signal*cos(2*pi*f_o) in frequency
my_fft2=fft(my2);
my_current2=abs(my_fft2)*1/mf_c;
mf=1:length(my_current2)

figure
plot(mf/length(my_current2)*mf_c,my_current2)

grid
legend('FFT')
axis([0 mf_c 0 max(my_current2)])



Il primo segnale raffigurato è un file wav. È il frutto di una manipolazione del segnale originale. Dopo una moltiplicazione per un coseno nel dominio del tempo è stato inserito un rumore cosinosuidale con uguale frequenza a quella di modulazione. Lo scopo è riottenere il segnale iniziale. Per farlo dovrei rimoltiplicarlo per lo stesso coseno della modulazione e filtrare fra -3.4 kHz e 3.4 kHz ricordando di eliminare la frequenza intorno allo zero.
Il problema è quando moltiplico il segnale per il coseno. Il risultato è... strano. Qualche buona anima può darmi una dritta :)?

[Russell]

credo che il tuo traslare lo spettro di 7500Hz non corrisponda alla frequenza di modulazione, e per questo il segnale dopo la tua manipolazione non torna ad essere quello originario a cavallo della frequenza 0Hz

Forse potrebbe aiutarti andare a disegnare lo spettro (FFT) specificando correttamente i valori dell'asse delle ascisse. Attualmente hai messo un contatore dei campioni, o qualcosa del genere.... ti suggerirei di ricavarti esattamente tale asse in Hz, a quel punto dovresti poter ragionare meglio sui numeri che stai utilizzando

[fairyvilje]

Al contrario, ho cercato di prestare attenzione e rappresentare l'asse in Hz. Se non fosse così vorrebbe dire che ho frainteso completamente il modo in cui si cambia scale D:
Il che è molto plausibile.

[Russell]

allora, dunque... scorrendo il tuo codice c'è questa riga che non mi piace

plot(mf/length(my_current)*32000,my_current)

mf/length(my_current) è un numero tra 0 e 1 e serve a scorrere l'asse delle ascisse da sinistra (0) a destra (1)

tu moltiplichi questo fattore per 32000... ottenendo frequenze da 0 a 32000
ma questo è un errore
infatti dopo la FFT hai un array con le frequenze da -fn a 0 a +fn (dove con fn ho indicato la frequenza di Nyquist = frequenza di campionamento /2)

segue che io scriverei piuttosto qualcosa tipo

plot(((mf/length(my_current))-0.5)*32000,my_current)

tanto per avere frequenze da -16000 a + 16000
(ecco che lo spettro diventa simmetrico rispetto all'origine)

non ho testato pero'... ragionaci da solo un attimo


EDIT. Adesso capisco da dove hai pescato quel 7500 Hz, hai osservato lo spettro FFT e hai guardato dove era il picco, stabilendone la frequenza ad occhio.... peccato che l'asse delle ascisse è errato...
Ad occhio, prendendo due misure spartane, la frequenza di modulazione è su 8500Hz

[fairyvilje]

Ora leggo il tuo post, volevo solo segnalare che se usassi 7499 Hz uscirebbe qualcosa mooolto più vicino a ciò che vorrei... E questa cosa mi spiazza...
Mi riguardo le tue note :)

[fairyvilje]

ma questo è un errore
infatti dopo la FFT hai un array con le frequenze da -fn a 0 a +fn (dove con fn ho indicato la frequenza di Nyquist = frequenza di campionamento /2)

Guarda anche secondo me, solo che in laboratorio ci hanno fatto rappresentare i moduli delle fft in questo modo. Quando ho chiesto perché... lasciamo stare

[fairyvilje]

Ecco il riferimento.

[fairyvilje]

infatti dopo la FFT hai un array con le frequenze da -fn a 0 a +fn (dove con fn ho indicato la frequenza di Nyquist = frequenza di campionamento /2)

Sei sicuro?
http://www.phys.nsu.ru/cherk/fft.pdf
Alla fine viene introdotta una funzione che fa quanto dici, fftshift.

[fairyvilje]

Sono abbastanza certo di aver trovato la soluzione del mio problema. Si chiama instabilità numerica D:.
Sono anche riuscito a sentire il messaggio segreto sfruttando il fatto che una scheda audio funge naturalmente da filtro passabanda .

Codice: Seleziona tutto

clear all
close all
clc

mf_c=32*10^3;

my=wavread("./Dati/segnaleDistorto_1.wav");
mt=1:length(my)

% Draw the signal in time
figure
plot(mt*1/mf_c,my)

grid
legend('Raw Wav')
axis([0 length(mt)*1/mf_c -1 1])


% Draw the signal in frequency
my_fft=fft(my);
my_current=fftshift(abs(my_fft))*1/mf_c;
mf=(-length(my_current)/2):(length(my_current)/2-1)

figure
plot((mf/length(my_current))*mf_c,my_current)

grid
legend('FFT')
axis([-mf_c/2 mf_c/2 0 max(my_current)])

[a,max_my_fft]=max(my_current);
max_my_fft=(max_my_fft-length(my_current)/2)/length(my_current)*mf_c


% Draw the signal in time
figure

for c=1:length(my)
  my2(c)=2*my(c)*cos(2*pi*1/mf_c*c*max_my_fft);
  c++
endfor

plot(mt*1/mf_c,my2);

grid
legend('Raw Wav')
axis([0 length(mt)*1/mf_c -1 1])


% Draw the signal in frequency
my_fft2=fft(my2);
my_current2=fftshift(abs(my_fft2))*1/mf_c;
mf=(-length(my_current2)/2):(length(my_current2)/2-1)

figure
plot(mf/length(my_current2)*mf_c,my_current2)

grid
legend('FFT')
axis([-mf_c/2 mf_c/2 0 max(my_current)])

wavwrite(my2,mf_c,"./Dati/segnaleDistorto_1a.wav");


Non riesco a convincermi degi 8500 Hz. Viene comunque 7500.

[Russell]

si, scusa ma mi sono un po' ingannato nel guardare il tuo spettro... e un po' dalla stanchezza di ieri sera
allora

le frequenze dopo la FFT in effetti partono da 0 (come hai concluso), solo che andrebbe sostanzialmente tagliata via la seconda parte del grafico perche' sarebbero le frequenze negative, ribaltate.... a tagliare via la seconda metà di un array si fa' presto in effetti con Matlab/Octave....altrimenti fare quei plot disegnando uno spettro fino a 32000Hz sarebbe un po' bruttino (avresti disegnato anche la parte in aliasing)
in pratica le frequenze dell'array che vedi dovrebbero essere: 0.... fn -fn .... 0

la fftshift è una funzione carina che si applica al segnale ottenuto dalla FFT: in pratica evita di incasinarsi con l'asse delle frequenze (vedi sopra) e all'atto pratico ribalta lo spettro stesso portandolo da "0.... fn -fn .... 0" a "-fn.... 0 .... fn", a quel punto diventa banale creare un array delle frequenze per fare correttamente i plot (sarebbe la formula che ti ho scritto io nel post precedente in pratica)

questo pezzo di codice dal tuo pdf dovrebbe essere tenuto a promemoria sempre:

X
=
abs(fft(x1,N));
X
=
fftshift(X);
F
=
[-N/2:N/2-1]/N;
plot(F,X)

detto questo quindi, hai ragione, la modulazione non dovrebbe essere avvenuta a 8500 Hz, ma 7500Hz come hai stimato

evidentemente, si, ci sara' qualche problema numerico da tenere di conto... o forse il modulatore a 7500Hz non era tanto preciso, come pure magari non era tanto precisa la frequenza di 32000Hz... bisogna in tal caso adattarsi, lo scostamento in frequenza in effetti è dell'ordine dell' 1 per-mille