MATLAB - Fast Fourier Transform FFT

[edomar]

Salve a tutti, vi espongo il mio problema nella speranza che possiate aiutarmi.

Dispongo di un vettore corrente ricavato rispetto a theta = omega*t (fondamentale 50 Hz):

Devo effettuare l'analisi in frequenza per poi ricavare l'armonica fondamentale (per calcolare il THD).

Ho seguito questo procedimento:

Codice: Seleziona tutto

% analisi armonica corrente i1

f_fond = 50;
omega = 2*pi*f_fond;
t = theta/omega;

TS = t(2) - t(1);
FS = 1/TS;

Y = fft(i1);

f = FS*(0:length(i1)-1)/length(i1);

stem(f,abs(Y)/length(i1)*2)
xlim([0 500])

Ottenendo questo spettro in ampiezza:

ZOOM

che è palesemente sbagliato. Qualcuno può aiutarmi? è la prima volta che mi cimento in un'analisi di questo tipo e sono sicuro di star facendo tante stupidaggini

[edomar]

EDIT: ho cambiato un po' il codice, rileggendo un po' di teoria sulla FFT... La situazione è migliorata, ma non capisco perché non mi rappresenti la fondamentale e le sue multiple, sembra tutto "shiftato".. Posto le modifiche:

Codice: Seleziona tutto

f_fond = 50;
omega = 2*pi*f_fond;
t = theta/omega;

TS = t(2) - t(1);
FS = 1/TS;

n = 2^nextpow2(length(i1))
Y = fft(i1,n);

f = FS*(0:n-1)/n;

stem(f,abs(Y)/n*2)
xlim([0 500])




Inoltre mi accorgo che variando il numero di elementi di theta (e quindi di t) il risultato cambia e anche di molto..

[gvee]

Perché posti un codice incompleto ?
Esempi per calcolare la DFT con MATLAB li puoi trovare in questo libro, p.e. al capitlo DFT Applications.

[edomar]

Ci sono 600 righe di codice prima... Posso postare "quello che serve"... Il problema è che non riesco a "estrarre" dal workspace il vettore i1.. Sai come si fa?

In rete ci sono tanti esempi, è vero... Ma sono tutti su combinazioni di sinusoidi, partendo dalla definizione di frequenza di campionamento, vettore del tempo, ecc ecc...
In questo caso ho un codice in cui dopo mille peripezie trovo questa corrente e devo studiarne gli spettri... Quindi l'intervallo di tempo è "vincolato"

[Maranza]

Help FFTshift!

[edomar]

Viene fuori una cosa improbabile... posto

Codice: Seleziona tutto

% analisi armonica corrente i1

omega = 2*pi*50;

t = theta/omega;

TS = t(2) - t(1);

FS = 1/TS;

n = 2^nextpow2(length(i1))

Y = fftshift(i1);

f = FS*(0:length(Y)-1)/length(Y);

stem(f,abs(Y)/length(Y)*2)
xlim([0 500])




Non capisco...
Sai come scrivere qui sul forum un vettore che ho salvato solo sul workspace? In questo modo posso postare theta e i1..

[gvee]

Sai come scrivere qui sul forum un vettore che ho salvato solo sul workspace? In questo modo posso postare theta e i1..

Prova a salvare i vettori in un file .mat ed allegarli al forum.
Oppure i un file .dat e ne copi qui il contenuto tra i tag code.

Comunque se conosci il vettore di tempo dovresti conoscere direttamente la frequenza di campionamento.

[edomar]

Esatto, è definita come TS nello script che ho postato.. Comunque il file .mat non è allegabile, il .dat mi dà fuori caratteri aramaici, tutto tranne il vettore.. non capisco il perché

[xyz]

Se definisci TS in questo modo non ha senso:

Codice: Seleziona tutto

...
TS = t(2) - t(1);
....

Il tempo di campionamento non può essere uguale alla differenza dei primi 2 campioni.

Comunque non credo che serva allegare un inutile matrice, basta che parti con un segnale generato nel codice facile come una armonica ad una determinata frequenza e poi vedi cosa ottieni facendo FFT. Se il risultato ti sembra corretto al posto della singola armonica gli dai il tuo segnale e vedi cosa ottieni.

[edomar]

Il vettore theta è equispaziato, per questo l'ho definito in questo modo..

hai ragione comunque, e il risultato è pure totalmente sbagliato:

Codice: Seleziona tutto

% analisi armonica corrente i1
theta = 0:0.001:2*pi;

omega = 2*pi*50;

t = theta/omega;

F = 10*sin(3*omega*t);

TS = t(2) - t(1);

FS = 1/TS;

n = 2^nextpow2(length(F));

Y = fft(F,n);

f = FS*(0:n-1)/n;

stem(f,abs(Y)/n*2)
xlim([0 500])