MATLAB FFT XAXIS límites que ensucian y fftshift

Question

Esta es la primera vez que estoy usando la función FFT y yo estoy tratando de trazar el espectro de frecuencia de una función coseno simple:

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

La frecuencia de muestreo es 220500. Estoy trazando un segundo de la función f.

Aquí es mi intento:

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]); 

¿Por qué recibo resultados extraños cuando aumente la frecuencia de 900 Hz? Estos resultados extraños pueden solucionarse aumentando los límites del eje x, por ejemplo, de 500 Hz a 1000 Hz. Además, ¿es este el enfoque correcto? Noté que muchas otras personas no usaron fftshift(X) (pero creo que solo hicieron un análisis de espectro de un solo lado).

Gracias.

Answer 1

Aquí está mi respuesta tal como prometió.

La primera o sus preguntas relacionadas con por qué "obtiene resultados impares cuando aumenta la frecuencia a 900 Hz" están relacionadas con la funcionalidad de cambio de escala de Matlab descrita por @Castilho. Cuando cambia el rango del eje x, Matlab intentará ser útil y reescalar el eje y. Si los picos se encuentran fuera de su rango especificado, matlab hará un acercamiento de los pequeños errores numéricos generados en el proceso. Puedes remediar esto con el comando 'ylim' si te molesta.

Sin embargo, su segunda pregunta, más abierta, "¿es este el enfoque correcto?" requiere una discusión más profunda. Permítame decirle cómo haría para hacer una solución más flexible para lograr su objetivo de trazar una onda de coseno.

se comienza con lo siguiente:

time = 1; 
freq = 220500; 

Esto plantea una alarma en mi cabeza inmediatamente. Mirando el resto de la publicación, parece que le interesan las frecuencias en el rango sub-kHz. Si ese es el caso, entonces esta tasa de muestreo es excesiva ya que el límite de Nyquist (sr/2) para esta tasa es superior a 100 kHz.Supongo que querías usar la frecuencia de muestreo de audio común de 22050 Hz (¿pero podría estar equivocado aquí?

De cualquier manera, su análisis funciona numéricamente OK al final. Sin embargo, no se está ayudando a sí mismo a comprender cómo la FFT se puede utilizar de manera más efectiva para el análisis en situaciones del mundo real.

Permíteme publicar cómo lo haría. La siguiente secuencia de comandos hace casi exactamente lo que hace la secuencia de comandos, pero abre un potencial sobre el que podemos construir. .

%// 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') 

Se calcula un eje de tiempo y se calcula el número de puntos de FFT a partir de la longitud del eje de tiempo. Esto es muy extraño. El problema con este enfoque es que la resolución de frecuencia de fft cambia a medida que cambia la duración de su señal de entrada, porque N depende de su variable de "tiempo". El comando matlab fft usará un tamaño FFT que coincida con el tamaño de la señal de entrada.

En mi ejemplo, calculo el eje de frecuencia directamente desde el NFFT. Esto es algo irrelevante en el contexto del ejemplo anterior, ya que configuro el NFFT para que sea igual al número de muestras en la señal. Sin embargo, usar este formato ayuda a desmitificar tu pensamiento y se vuelve muy importante en mi próximo ejemplo.

** NOTA LATERAL: Utiliza real (F) en su ejemplo. A menos que tenga una muy buena razón para extraer solo la parte real del resultado de FFT, entonces es mucho más común extraer la magnitud de la FFT usando abs (F). Este es el equivalente de sqrt (real (F).^2 + imag (F).^2). **

La mayoría de las veces querrá usar un NFFT más corto. Esto podría deberse a que quizás está ejecutando el análisis en un sistema de tiempo real, o porque desea promediar el resultado de muchas FFT para tener una idea del espectro promedio de una señal variable en el tiempo, o porque desea comparar los espectros de señales que tienen diferente duración sin desperdiciar información. Simplemente usando el comando fft con un valor de NFFT <, la cantidad de elementos en su señal dará como resultado un fft calculado a partir de los últimos puntos NFFT de la señal. Esto es un poco derrochador

El siguiente ejemplo es mucho más relevante para aplicaciones útiles. Se muestra cómo se divide una señal en bloques y luego procesar cada bloque y promediar el resultado:

%//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') 

utilizo una ventana de Hamming en el ejemplo anterior. La ventana que elija debe adaptarse a la aplicación http://en.wikipedia.org/wiki/Window_function

La cantidad de solapamiento que elija dependerá en cierto modo del tipo de ventana que utilice. En el ejemplo anterior, la ventana de Hamming pondera las muestras en cada memoria intermedia hacia cero desde el centro de cada cuadro. Para usar toda la información en la señal de entrada, es importante usar un poco de superposición. Sin embargo, si solo usa una ventana rectangular simple, la superposición no tiene sentido ya que todas las muestras tienen el mismo peso. Cuanta más superposición utilice, más procesamiento se requiere para calcular el espectro medio.

Espero que esto ayude a su comprensión.

Answer 2

Su resultado es perfecto. El cálculo del eje de frecuencia también es correcto. El problema radica en la escala del eje y. Cuando utiliza la función xlims, matlab recalcula automáticamente la escala y para que pueda ver datos "significativos". Cuando los picos del coseno se encuentran fuera del límite que usted eligió (cuando f> 500Hz), no hay picos que mostrar, por lo que la escala se calcula en función de algún pequeño ruido (aquí en mi computadora, con matlab 2011a, la escala y fue 10 -dieciséis).

Cambiar el límite es de hecho el enfoque correcto, porque si no lo cambia no puede ver los picos en el espectro de frecuencia.

Una cosa que noté, sin embargo. ¿Hay alguna razón para que traces la parte real de la transformación? Por lo general, es abs(F) que se traza, y no la parte real.

editar: En realidad, su eje de frecuencia es correcto porque df, en este caso, es 1. La línea de faxes es correcta, pero el cálculo de df no es correcto.

La FFT calcula N puntos de -Fs/2 a Fs/2. Por lo tanto, N puntos sobre un rango de F rinde un df de Fs/N. Como N/tiempo = Fs => tiempo = N/Fs. Sustituyendo eso en la expresión de df que utilizó: your_df = Fs/N * (N/Fs) = (Fs/N)^2. Como Fs/N = 1, el resultado final fue el correcto: P