PROGRAMA BASIC PARA CALCULO DE LA FFT
Llevando a cabo las indicaciones del diagrama de flujo, se llega al siguiente programa BASIC que se combina con Mathcad para evaluación.
Se va a trabajar con el ejemplo de la exponencial decreciente a la que se le hallará la Transformada de Fourier.
Archivo ASCII que contiene la parte real de los datos
Archivo ASCII que contiene la parte imaginaria de los datos
Una vez que los datos se han incorporado en los archivos ent_real.prn y ent_imag.prn se ejecuta el siguiente programa BASIC:
5 pi = 4 * ATN(1): GOSUB 1000 'Pantalla
10 INPUT "Número de Muestras (N): "; N
15 INPUT "Periodo de muestreo (T): "; T
20 nu = INT(LOG(N) / LOG(2) + .0001)
30 DIM xr(N + 1), xi(N + 1), yr(N + 1), yi(N + 1)
40 ON a GOTO 2000, 3000, 4000, 5000
1000 CLS
1010 PRINT "Elija: 1 - Para fft"
1011 PRINT " 2 - Para ifft"
1012 PRINT " 3 - Para Convolución"
1013 PRINT " 4 - Para Correlación"
1020 INPUT a
1030 RETURN
2000 REM Rutina fft
2010 GOSUB 6000 'Ingreso de valores
2020 GOSUB 7000 'Núcleo fft
2030 GOTO 8000 'Graficación
3000 REM Rutina ifft
3010 GOSUB 6000
3020 FOR k = 0 TO N - 1: xi(k) = -xi(k): NEXT k
3030 GOSUB 7000
3040 FOR k = 0 TO N - 1: xr(k) = xr(k) / N: xi(k) = xi(k) / N: NEXT k
3050 GOTO 8000
4000 REM Rutina convolución
4010 GOSUB 6000 'Ingreso de valores de la primera Función
4020 GOSUB 7000 'Núcleo fft
4030 FOR k = 0 TO N - 1: yr(k) = xr(k): yi(k) = xi(k): NEXT k
4045 FOR k = 0 TO N - 1: xr(k) = 0: xi(k) = 0: NEXT k
4047 GOSUB 6000 'Ingreso de valores de la segunda Función
4050 GOSUB 7000 'Núcleo fft
4060 FOR k = 0 TO N - 1
4062 c1 = xr(k): xr(k) = xr(k) * yr(k) - xi(k) * yi(k)
4064 xi(k) = -c1 * yi(k) - yr(k) * xi(k)
4066 NEXT k
4070 GOSUB 7000
4080 FOR k = 0 TO N - 1: xr(k) = xr(k) / N: xi(k) = xi(k) / N: NEXT k
4090 GOTO 8000
5000 REM Rutina correlación
5010 GOSUB 6000 'Ingreso de valores de la primera Función
5020 GOSUB 7000 ' Núcleo fft
5030 FOR k = 0 TO N - 1: yr(k) = xr(k): yi(k) = xi(k): NEXT k
5040 FOR k = 0 TO N - 1: xr(k) = 0: xi(k) = 0: NEXT k
5050 GOSUB 6000 'Ingreso de valores de la segunda función
5060 GOSUB 7000 'Núcleo fft
4075 FOR k = 0 TO N - 1: yi(k) = -xi(k): NEXT k
5080 FOR k = 0 TO N - 1
5082 c1 = xr(k): xr(k) = xr(k) * yr(k) - xi(k) * yi(k)
5084 xi(k) = -c1 * yi(k) - yr(k) * xi(k)
5086 NEXT k
5110 GOSUB 7000 'Halla ifft
5120 FOR k = 0 TO N - 1: xr(k) = xr(k) / N: xi(k) = xi(k) / N: NEXT k
5130 GOTO 8000 'Rutina Graficación
6000 REM Subrutina de ingreso de valores reales
6004 OPEN "c:\winmcad\ent_real.prn" FOR INPUT AS #1
6005 FOR k = 0 TO N - 1: INPUT #1, xr(k): NEXT k
6008 CLOSE #1
6010 REM Subrutina de ingreso de valores imaginarios
6020 OPEN "c:\winmcad\ent_imag.prn" FOR INPUT AS #1
6030 FOR k = 0 TO N - 1: INPUT #1, xi(k): NEXT k
6040 CLOSE #1
'6010 FOR k = 0 TO N - 1
'6020 PRINT "valor real "; k; " de la entrada": INPUT xr(k)
'6030 NEXT k
'6040 FOR k = 0 TO N - 1
'6050 PRINT "valor imaginario "; k; " de la entrada": INPUT xi(k)
'6060 NEXT k
6070 RETURN
7000 REM N£cleo fft
7010 L = 1: n2 = N / 2: n1 = nu - 1: k = 0
7020 IF L > nu THEN 7140
7030 i = 1
7040 M = INT(k / 2 ^ n1)
7050 i2 = 0: ib = 0: z = M: GOSUB 7200: P = ib
7060 ar = 2 * pi * P / N: c = COS(ar): s = SIN(ar): ag = k + n2
7070 tr = xr(ag) * c + xi(ag) * s: ti = xi(ag) * c - xr(ag) * s
7080 xr(ag) = xr(k) - tr: xi(ag) = xi(k) - ti
7090 xr(k) = xr(k) + tr: xi(k) = xi(k) + ti: k = k + 1
7100 IF i = n2 THEN 7110 ELSE i = i + 1: GOTO 7040
7110 k = k + n2
7120 IF k < N - 1 THEN 7030
7130 k = 0: n1 = n1 - 1: n2 = n2 / 2: L = L + 1: GOTO 7020
7140 i2 = 0: ib = 0: z = k: GOSUB 7200: i = ib
7150 IF i <= k THEN 7180
7160 t3 = xr(k): t4 = xi(k): xr(k) = xr(i): xi(k) = xi(i)
7170 xr(i) = t3: xi(i) = t4
7180 IF k = N - 1 THEN 7190 ELSE k = k + 1: GOTO 7140
7190 RETURN
7200 IF i2 < nu THEN 7210 ELSE RETURN
7210 j2 = INT(z / 2): ib = 2 * ib + z - 2 * j2: z = j2: i2 = i2 + 1: GOTO 7200
8000 REM Rutina Graficación
8002 OPEN "c:\winmcad\salida.prn" FOR OUTPUT AS #1
8004 FOR k = 0 TO N - 1: PRINT #1, xr(k) * T, xi(k) * T: NEXT k
8006 CLOSE #1: END
'8003 FOR k = 0 TO N - 1
'8005 yr(k) = SQR(xr(k) ^ 2 + xi(k) ^ 2)
'8006 PRINT USING "##.##"; xr(k); : PRINT " "; : PRINT USING "##.##"; xi(k)
'8007 NEXT k: END
8008 SCREEN 2: M = xr(0)
8010 FOR k = 0 TO N - 1
8020 IF M >= xr(k) THEN 8040
8030 M = xr(k)
8040 NEXT k
8050 q = 90 / M
8060 FOR k = 0 TO 255: PSET (k, 90): NEXT k
8070 FOR k = 0 TO N - 1: PSET (k * 245 / N, q * xr(k)): NEXT k
8105 j$ = INKEY$
8110 IF j$ = "" THEN 8105
8120 FOR k = 0 TO N - 1: xr(k) = yr(k): NEXT k: GOTO 8008
Una vez corrido el programa, devuelve el archivo salida.prn una matriz de dos columnas, donde la primera corresponde a la parte real de la transformada y la segunda a la parte imaginaria.
Si se observa la parte real:
Y la parte imaginaria
Lo que coincide plenamente con el cálculo directo de la Transformada de Fourier.
TRANSFORMADA DE FOURIER CON MATHCAD
fft(v)
Retorna la TF de un vector. El resultado es un vector de 1+2m-1 elementos cuyo n-ésimo elemento está dado por
Argumentos:
- v debe tener 2m elementos (m > 2).
- Todos los elementos en v son reales.
Para vectores con valores complejos o con cualquier número de elementos, use cfft en lugar de fft.
Para TF bi-dimensionales , use cfft.
_____________________________________________________________________
ifft(v)
Retorna la IFT de un vector. El resultado es un vector de 1+2m-1 elementos cuyo n-ésimo elemento está dado por
_____________________________________________________________________
cfft(A)
Retorna la FT de un vector o matriz. El resultado tiene el mismo número de filas y columnas que A. Si A es un vector, el resultado está dado por:
Argumentos:
- A puede ser un vector o una matíz.
Si...
- A es un vector tiene 2m elementos, y..
- Todos los elementos en el vector son reales
...use fft en lugar de cfft. Para tales vectores, la segunda mitad del espectro es la imagen del espejo de la primera mitad y no necesitan ser calculados.
Para argumentos matriciales, cfft retorna la FT bidimensional:
_____________________________________________________________________
icfft(B)
Retorna la IFT correspondiente a la fft.
Argumentos:
B puede ser una matriz o un vector.
_____________________________________________________________________
FFT(v)
Retorna la FT de un vector. El resultado es un vector de 1+2m-1 elementos cuyo n-ésimo elemento está dado por
Argumentos:
- v debe tener 2m elementos (m > 2).
- Todos los elementos en v son reales.
Para vectores con valores complejos o con cualquier número de elementos, use CFFT en lugar de FFT.
Para FT bi-dimensionales , use CFFT.
_____________________________________________________________________
IFFT(u)
Retorna la IFT de un vector correspondiente a FFT. El resultado es un vector de 2m elementos cuyo n-ésimo elemento está dado por
Argumentos:
- u es un vector con 1+2m-1 elementos (m > 2).
- Todos los elementos en u son reales.
Se debe usar esta función sólo cuando los datos sean puramente reales. Si los datos son complejos se debe usar ICFFT.
CFFT(A)
Retorna la FT de un vector o matriz. El resultado tiene el mismo número de filas y columnas que A. Si A es un vector, el resultado está dado por:
Argumentos:
- A puede ser un vector o una matriz.
Si...
- Todos los elementos en el vector son reales
...use FFT en lugar de CFFT. Para tales vectores, la segunda mitad del espectro es la imagen del espejo de la primera mitad y no necesitan ser calculados.
Para argumentos matriciales, CFFT retorna la FT bidimensional:
_____________________________________________________________________
ICFFT(B)
Retorna la IFT correspondiente a la CFFT.
Argumentos:
B puede ser una matriz o un vector.
_____________________________________________________________________
TOPICOS ADICIONALES
Frecuencias correspondientes a los coeficientes:
Tanto fft como cfft retornan vectores cuyos elementos son las amplitudes complejas de las distintas frecuencias comprendidas en la señal original.
Para recuperar dichas frecuencias, se debe conocer:
- La frecuencia de muestreo de la señal original.
- El número de muestras en la señal original.
Dados estos parámetros, la frecuencia asociada con el n-ésimo elemento está dada por:
Donde: fs es la frecuencia de muestreo y N es el número de muestras.
Aliasing:
Debido a que la fft es una aproximación discreta a la transformada continua, se pueden encontrar armónicas espúreas. Esto se llama Aliasing. Para evitar el aliasing:
- Asegurarse que la señal tiene ancho de banda finito
Espectro de Potencia y de Fase:
El Espectro de Potencia en el dominio de frecuencia está dado por:
El Espectro de Fase en el dominio de frecuencia está dado por:
Las funciones fft y cfft retornan las partes reales e imaginarias de la Transformada de Fourier. Se puede recuperar el espectro de potencia y el de fase a partir de estas usando el operador vectorización:
fft/ifft contra cfft/cifft:
Usar el par fft/ifft si ambas de las siguientes sentencias son ciertas:
- Los datos en el dominio del tiempo son reales.
- El vector de datos tiene 2m elementos.(m>2).
Usar el par cfft/icfft si cualquiera de las siguientes sentencias son ciertas:
- Los datos en el dominio dle tiempo son complejos.
- El vector de datos no tiene 2m elementos.
- Los datos son una matriz, no un vector.
Para datos reales en el dominio del tiempo, la FT tiene simetría conjugada. Debido a esto, fft caerá en una segunda mitad redundante. Este es el porqué el vector retornado por fft es de la mitad de tamaño del vector original.
La función ifft reconstruye esta mitad redundante. Por esta razón, se debería usar fft e ifft juntas. El par cfft/icfft no hace tales suposiciones acerca de simetría. Aunque estas funciones no son tan eficientes para datos reales, no obstante pueden ser usadas.
OPERACIONES POSIBLES DE REALIZAR CON EL ALGORITMO FFT
Ejemplo de Transformada de Fourier de una señal
Una vez corrido el programa, devuelve el archivo sal_fft.prn una matriz de dos
columnas, donde la primera corresponde a la parte real de la transformada y la segunda a la parte imaginaria.
Si se observa la parte real:
Lo que coincide plenamente con el cálculo directo de la Transformada de Fourier:
Si se analiza la parte imaginaria:
Ejemplo de Transformada Inversa de Fourier de una señal
Una vez corrido el programa, devuelve el archivo sal_ifft.prn una matriz de dos columnas, donde la primera corresponde a la parte real de la transformada Inversa y la segunda a la parte imaginaria.
Si se observa la parte real:
Lo que coincide plenamente con el cálculo directo de la Transformada Inversa de Fourier:
Si se analiza la parte imaginaria:
Ejemplo de Convolución de dos señales
Una vez corrido el programa, devuelve el archivo sal_conv.prn una matriz de dos columnas, donde la primera corresponde a la parte real de la transformada y la segunda a la parte imaginaria.
Si se observa la parte real:
Lo que coincide plenamente con el cáculo directo de la Convolución:
Ejemplo de Correlación
Si se observa la parte real:
Lo que coincide plenamente con el cálculo directo de la Correlación:
Si se tienen dos secuencias reales g(k.T) y h(k.T), ambas de longitud N, se puede encontrar las FFT simultáneamente formando una secuencia compleja a(k,T) = g(k.T) + j h(k.T). Usando esta secuencia como entrada al programa BASIC de cálculo, se obtienen las correspondientes FFT, A(k) la cual es es G(k) + j H(k). Después, la secuencia a*(k.T) = g(k.T) - j h(k.T) se forma y su FFT se determina. Esto está dado por:
A*(N-k) = G(k) - j H(-k)
donde G(k) y H(k) se extraen de A(k) y A*(N-k) formando las siguientes ecuaciones:
Volver a la página Anterior.