rapida - transformada inversa de fourier matlab

Transformada analítica de Fourier vs FFT de funciones en Matlab

Question

He adaptado el código en Comparando FFT of Function to Analytical FT Solution en Matlab para esta pregunta. Estoy tratando de hacer FFT y comparar el resultado con expresiones analíticas en las tablas de Wikipedia .

Mi código es:

a = 1.223;
fs = 1e5; %sampling frequency
dt = 1/fs;
t = 0:dt:30-dt;     %time vector
L = length(t); % no. sample points
t = t - 0.5*max(t); %center around t=0

y = ; % original function in time
Y = dt*fftshift(abs(fft(y))); %numerical soln

freq = (-L/2:L/2-1)*fs/L; %freq vector
w = 2*pi*freq; % angular freq

F = ; %analytical solution

figure; subplot(1,2,1); hold on
plot(w,real(Y),''.'')
plot(w,real(F),''-'')
xlabel(''Frequency, w'')
title(''real'')
legend(''numerical'',''analytic'')
xlim([-5,5])
subplot(1,2,2); hold on;
plot(w,imag(Y),''.'')
plot(w,imag(F),''-'')
xlabel(''Frequency, w'')
title(''imag'')
legend(''numerical'',''analytic'')
xlim([-5,5])

Si estudio la función gaussiana y dejo

y = exp(-a*t.^2); % original function in time

F = exp(-w.^2/(4*a))*sqrt(pi/a); %analytical solution

en el código anterior, parece que hay un buen acuerdo cuando se trazan las partes real e imaginaria de la función:

Pero si estudio un exponencial en descomposición multiplicado por una función Heaviside:

H = @(x)1*(x>0); % Heaviside function
y = exp(-a*t).*H(t);

F = 1./(a+1j*w); %analytical solution

entonces

¿Por qué hay una discrepancia? Sospecho que está relacionado con la línea Y = pero no estoy seguro de por qué o cómo.

Editar: cambié el ifftshift a fftshift en Y = dt*fftshift(abs(fft(y))); . Luego también eliminé los abs . El segundo gráfico ahora se ve así:

¿Cuál es la razón matemática detrás del gráfico ''reflejado'' y cómo puedo eliminarlo?

Answer 1

Las tramas al final de la pregunta no se reflejan. Si traza los que usan líneas en lugar de puntos, verá que los resultados numéricos tienen frecuencias muy altas. El componente absoluto coincide, pero la fase no. Cuando esto sucede, es casi seguro que se trate de un cambio en el dominio del tiempo.

Y, de hecho, usted define la función del dominio del tiempo con el origen en el medio. La FFT espera que el origen esté en la primera muestra (más a la izquierda). Para esto sirve ifftshift :

Y = dt*fftshift(fft(ifftshift(y)));

ifftshift mueve el origen a la primera muestra, en preparación para la llamada fft , y fftshift mueve el origen del resultado a la mitad, para mostrarlo.

Editar

Tu t no tiene un 0:

>> t(L/2+(-1:2))
ans =
  -1.5000e-05  -5.0000e-06   5.0000e-06   1.5000e-05

La muestra en t(floor(L/2)+1) debe ser 0. Esa es la muestra que ifftshift desplazamiento a la izquierda se mueve a la muestra más a la izquierda. (Uso el floor allí en caso de que L tenga un tamaño extraño, no el caso aquí).

Para generar una t correcta, haga lo siguiente:

fs = 1e5; % sampling frequency
L = 30 * fs;
t = -floor(L/2):floor((L-1)/2);
t = t / fs;

Primero genero un eje t entero de la longitud correcta, con 0 en la ubicación correcta ( t(floor(L/2)+1)==0 ). Luego lo convierto a segundos dividiéndolo por la frecuencia de muestreo.

Con esta t , la Y como sugiero anteriormente, y el resto de su código tal como está, veo esto para el ejemplo gaussiano:

>> max(abs(F-Y))
ans =    4.5254e-16

Para la otra función veo diferencias mayores, en el orden de 6e-6. Esto se debe a la incapacidad de muestrear la función Heaviside. Necesita t = 0 en su función muestreada, pero H no tiene un valor en 0. Noté que el componente real tiene un desplazamiento de magnitud similar, que es causado por la muestra en t = 0.

Típicamente, la función Heaviside muestreada se establece en 0.5 para t = 0 . Si hago eso, el desplazamiento se elimina por completo, y la diferencia máxima para el componente real se reduce en 3 órdenes de magnitud (los errores más grandes ocurren para valores muy cercanos a 0, donde veo un patrón de zig-zag). Para el componente imaginario, el error máximo se reduce a 3e-6, aún bastante grande, y es máximo a altas frecuencias. Atribuyo estos errores a la diferencia entre las funciones Heaviside ideales y muestreadas.

Probablemente debería limitarse a las funciones de banda limitada (o casi limitadas como la gaussiana). Es posible que desee intentar reemplazar la función Heaviside con una función de error (integral de Gaussian) con una pequeña sigma ( sigma = 0.8 * fs es la sigma más pequeña que consideraría para un muestreo adecuado). Su transformada de Fourier es conocida .