MATLAB FFT XAxis limiti rovinare e fftshift

Question

Questa è la prima volta che sto utilizzando la funzione FFT e sto cercando di tracciare lo spettro di frequenza di una semplice funzione coseno:

f = cos (2 * pi * 300 * t)

La frequenza di campionamento è 220500. Sto tracciando un secondo della funzione f.

Ecco il mio tentativo:

time = 1; 
freq = 220500; 
t = 0 : 1/freq : 1 - 1/freq; 
N = length(t); 
df = freq/(N*time); 

F = fftshift(fft(cos(2*pi*300*t))/N); 
faxis = -N/2/time : df : (N/2-1)/time; 

plot(faxis, real(F)); 
grid on; 
xlim([-500, 500]); 

Perché ricevo strani risultati quando ho aumentare la frequenza a 900Hz? Questi strani risultati possono essere risolti aumentando i limiti dell'asse x da, ad esempio, da 500 Hz a 1000 Hz. Inoltre, questo è l'approccio corretto? Ho notato che molte altre persone non hanno usato fftshift(X) (ma penso che abbiano fatto solo un'analisi spettrale su un lato).

Grazie.

Answer 1

Ecco la mia risposta come promesso.

La prima o le tue domande relative al motivo per cui "ottieni risultati strani quando aumenti la frequenza a 900 Hz" è correlata alla funzionalità di riscalamento della trama di Matlab come descritto da @Castilho. Quando modifichi l'intervallo dell'asse x, Matlab proverà ad essere utile e ridimensionare l'asse y. Se i picchi si trovano al di fuori dell'intervallo specificato, MATLAB ingrandirà i piccoli errori numerici generati nel processo. Puoi rimediare a questo con il comando 'ylim' se ti dà fastidio.

Tuttavia, la seconda, più aperta domanda "è questo l'approccio corretto?" richiede una discussione più profonda. Consentitemi di dirvi come farei per creare una soluzione più flessibile per raggiungere il vostro obiettivo di tracciare un'onda cosina.

Si inizia con il seguente:

time = 1; 
freq = 220500; 

Ciò solleva immediatamente un allarme nella mia testa. Guardando il resto del post, sembra che ti interessi le frequenze nella gamma sub-kHz. In questo caso, questa frequenza di campionamento è eccessiva poiché il limite di Nyquist (sr/2) per questa frequenza è superiore a 100 kHz.Immagino che intendiate usare la frequenza di campionamento audio comune di 22050 Hz (ma potrei sbagliarvi qui)?

In entrambi i casi, l'analisi funziona numericamente OK alla fine. Tuttavia, non stai aiutando te stesso a capire come la FFT può essere utilizzata in modo più efficace per l'analisi in situazioni reali.

Consentitemi di postare come farei questo. Il seguente script fa esattamente ciò che fa lo script, ma apre alcune potenzialità su cui possiamo costruire. .

%// These are the user parameters 
durT = 1; 
fs = 22050; 
NFFT = durT*fs; 
sigFreq = 300; 

%//Calculate time axis 
dt = 1/fs; 
tAxis = 0:dt:(durT-dt); 

%//Calculate frequency axis 
df = fs/NFFT; 
fAxis = 0:df:(fs-df); 

%//Calculate time domain signal and convert to frequency domain 
x = cos( 2*pi*sigFreq*tAxis ); 
F = abs( fft(x, NFFT)/NFFT ); 

subplot(2,1,1); 
plot( fAxis, 2*F ) 
xlim([0 2*sigFreq]) 
title('single sided spectrum') 

subplot(2,1,2); 
plot( fAxis-fs/2, fftshift(F) ) 
xlim([-2*sigFreq 2*sigFreq]) 
title('whole fft-shifted spectrum') 

Si calcola un asse del tempo e si calcola il numero di punti FFT dalla lunghezza dell'asse del tempo. Questo è molto strano Il problema con questo approccio è che la risoluzione in frequenza del fft cambia al cambiare della durata del segnale in ingresso, perché N dipende dalla variabile "tempo". Il comando matlab fft utilizzerà una dimensione FFT che corrisponde alla dimensione del segnale di ingresso.

Nel mio esempio, calcolo l'asse della frequenza direttamente dall'NFFT. Questo è alquanto irrilevante nel contesto dell'esempio precedente, poiché imposto che l'NFFT sia uguale al numero di campioni nel segnale. Tuttavia, l'uso di questo formato aiuta a demistificare il tuo pensiero e diventa molto importante nel mio prossimo esempio.

** NOTA LATO: si utilizza reale (F) nell'esempio. A meno che tu non abbia una buona ragione per estrarre solo la parte reale del risultato FFT, allora è molto più comune estrarre la grandezza della FFT usando abs (F). Questo è l'equivalente di sqrt (real (F).^2 + imag (F).^2). **

La maggior parte delle volte si vorrà utilizzare un NFFT più corto. Ciò potrebbe essere dovuto al fatto che si sta eseguendo l'analisi in un sistema in tempo reale o perché si desidera calcolare la media del risultato di molti FFT per ottenere un'idea dello spettro medio per un segnale che varia nel tempo o perché si desidera confrontare gli spettri di segnali che hanno durata diversa senza sprecare informazioni. Utilizzando il comando fft con un valore di NFFT <, il numero di elementi nel segnale genererà un valore calcolato dagli ultimi punti NFFT del segnale. Questo è un po 'dispendioso.

Il seguente esempio è molto più rilevante per l'applicazione utile. Essa mostra come si potrebbe dividere un segnale in blocchi e poi elaborare ogni blocco e media il risultato:

%//These are the user parameters 
durT = 1; 
fs = 22050; 
NFFT = 2048; 
sigFreq = 300; 

%//Calculate time axis 
dt = 1/fs; 
tAxis = dt:dt:(durT-dt); 

%//Calculate frequency axis 
df = fs/NFFT; 
fAxis = 0:df:(fs-df); 

%//Calculate time domain signal 
x = cos( 2*pi*sigFreq*tAxis ); 

%//Buffer it and window 
win = hamming(NFFT);%//chose window type based on your application 
x = buffer(x, NFFT, NFFT/2); %// 50% overlap between frames in this instance 
x = x(:, 2:end-1); %//optional step to remove zero padded frames 
x = ( x' * diag(win) )'; %//efficiently window each frame using matrix algebra 

%// Calculate mean FFT 
F = abs( fft(x, NFFT)/sum(win) ); 
F = mean(F,2); 

subplot(2,1,1); 
plot( fAxis, 2*F ) 
xlim([0 2*sigFreq]) 
title('single sided spectrum') 

subplot(2,1,2); 
plot( fAxis-fs/2, fftshift(F) ) 
xlim([-2*sigFreq 2*sigFreq]) 
title('whole fft-shifted spectrum') 

Io uso una finestra di Hamming nell'esempio di cui sopra. La finestra scelta deve essere adatta all'applicazione http://en.wikipedia.org/wiki/Window_function

La quantità di sovrapposizione scelta dipende in qualche modo dal tipo di finestra che si utilizza. Nell'esempio sopra, la finestra Hamming pesa i campioni in ciascun buffer verso zero lontano dal centro di ciascun fotogramma. Per utilizzare tutte le informazioni nel segnale di ingresso, è importante utilizzare alcune sovrapposizioni. Tuttavia, se si utilizza semplicemente una finestra rettangolare semplice, la sovrapposizione diventa inutile poiché tutti i campioni sono ugualmente pesati. Più si sovrappongono, più l'elaborazione è necessaria per calcolare lo spettro medio.

Spero che questo aiuti la vostra comprensione.

Answer 2

Il risultato è perfetto. Anche il calcolo dell'asse di frequenza è corretto. Il problema si trova sulla scala dell'asse y. Quando si utilizza la funzione xlims, matlab ricalcola automaticamente la scala y in modo da poter visualizzare dati "significativi". Quando i picchi del coseno giacciono al di fuori del limite che hai scelto (quando f> 500Hz), non ci sono picchi da mostrare, quindi la scala è calcolata sulla base di qualche piccolo rumore (qui al mio computer, con matlab 2011a, la scala y era 10 -16).

Cambiare il limite è effettivamente l'approccio corretto, perché se non lo si modifica non è possibile vedere i picchi sullo spettro di frequenza.

Una cosa che ho notato, tuttavia. C'è una ragione per te per tracciare la parte reale della trasformazione? Di solito, è il abs(F) che viene tracciato, e non la parte reale.

modifica: In realtà, l'asse delle frequenze è corretto solo perché df, in questo caso, è 1. La linea fax è corretta, ma il calcolo df non lo è.

La FFT calcola N punti da -Fs/2 a Fs/2. Quindi N punti su un intervallo di Fs produce un df di Fs/N. As N/time = Fs => time = N/Fs. Sostituendolo con l'espressione di df che hai usato: your_df = Fs/N * (N/Fs) = (Fs/N)^2. Come Fs/N = 1, il risultato finale era corretto: P