Source: https://www.scribd.com/document/75445987/AnalisisEspectralMedianteUsoFFT
Timestamp: 2018-09-24 23:02:56+00:00

Document:
Uploaded by lucho999
laboratorio 3 señales
Informe-Diseno-SisAdqu
Pfc 3181
Informatica II (Plantel 17) (2)
Comandos Básicos Para el Manejo de Señales en Matlab.docx
Practica Capetito
Guia Integrada de Actividades Academicas Ver Diciembre 12 2015
Mendez Conny Numero Log i A
Análisis espectral mediante el uso de la FFT
Mario Estévez Báez1 Andrés Machado García2 José M. Estévez Carrera3
Material publicado originalmente en formato html en: librosabiertos:analisis_espectral_mediante_el_uso_de_la_fft. InfoWiki. January 8, 2008, 09:08 CST. Available at: http://infomed20.sld.cu/wiki/doku.php?id=librosabiertos:analisis_espectral_mediante_el_uso_de_la_fft&rev=1199801296. Accessed January 8, 2008
El análisis espectral constituye un método de análisis poderoso para poner en evidencia periodicidades ocultas en una serie temporal, tal como lo es, por ejemplo, el electroencefalograma. El espectro de potencia refleja la energía de cada uno de los componentes de frecuencia del proceso estudiado. Permite distinguir los componentes espectrales y cuantificarlos. Este método de análisis en Neurofisiología Clínica, cuenta ya con más de 60 años de aplicación, existiendo equipos comerciales que permiten su empleo para el estudio de la actividad bioeléctrica del cerebro. Éste método de análisis es conocido y se maneja de forma diaria por el especialista promedio. En el caso del electroencefalograma, el espectro de interés se encuentra entre 0.5 y 50 hertzios, mientras que para aquellas aplicaciones que muestren las fluctuaciones de variables como el neumograma, la presión arterial y el ritmo cardíaco, por poner solo algunos ejemplos, se requiere analizar frecuencias menores de 0.5 hertzios. Teniendo ello en cuenta, resulta conveniente revisar algunos aspectos básicos del método en cuestión, y en donde sea necesario se realizarán referencias específicas, vinculadas con el estudio de tales frecuencias, que por su carácter han sido denominadas extra-lentas y que están vinculadas con manifestaciones de fenómenos vitales del organismo del hombre y los animales (Estévez Báez M.., 1982).
Doctor en Medicina, Especialista de Fisiología de Segundo Grado, Investigador Titular, Profesor Consultante, Doctor en Ciencias Médicas, Académico Titular AIA, Instituto de Endocrinología y Enfermedades Metabólicas MINSAP. 2 Licenciado en Cibernética-Matemática, Profesor Auxiliar, Maestro en Ciencias de la Computación Facultad de Biología, Universidad de La Habana, MES. 3 Licenciado en Informática, Instituto Superior de Medicina Militar “Dr. Luis Díaz Soto”
La secuencia temporal x(n), que es periódica, con período “T”, puede ser representada, como hemos visto en otro acápite anterior, por la siguiente expresión:
x(n) = x( n + rT ) ,
donde "r" es un valor entero cualquiera. Esta secuencia puede ser representada por una serie de Fourier, que corresponda a una suma de secuencias exponenciales complejas relacionadas, con frecuencias que sean múltiplos de la frecuencia fundamental (2π/T), asociada con la secuencia periódica x(n). Las exponenciales periódicas complejas son de la siguiente forma: e i ( 2 π T ) kn Donde "i" es la unidad imaginaria y "k" es un número entero. La representación de la serie de Fourier de la secuencia temporal x(n) toma la siguiente forma:
1 T −1 ∑0 x(k ) e i ( 2π T ) kn . T k=
En esta expresión, queda bien definido que se trata de una secuencia de una señal de valores discretos y que por lo tanto, requiere solamente de “T” exponenciales complejas relacionadas. Los coeficientes de la serie discreta de Fourier, pueden ser obtenidos de la secuencia temporal x(n), por la relación que se muestra en la siguiente expresión:
x(k ) = ∑ x(n) e −i ( 2 π T ) kn ;
En esta expresión se puede generalizar que
x(k ) = x (T + k ) ,
para cualquier entero “k”.
Las expresiones "a" y "b" constituyen una pareja (par) de análisis-síntesis, conociéndose como las expresiones de las Series Discretas de Fourier, para la representación de una secuencia periódica. Si se sustituye, para facilitar la notación, al término exponencial complejo de la siguiente manera:
WT = e − i ( 2 π T ) ,
sustracciones y adiciones. y tomando en cuenta que la serie temporal tiene que ser compleja. para secuencias temporales de tiempo discreto. n=0 T −1 En la anterior exposición. para el cálculo de una serie temporal periódica de tiempo finito: X (k ) = T −1 n=0 ∑ X ( n) e T −1 n=0 −i ( 2π / T ) kn que se puede simplificar como fue descrito anteriormente de la siguiente manera: X (k ) = ∑ X (n)WT kn .. para una secuencia temporal discreta. en forma de exponenciales complejas. análisis : x(k ) = ∑ x(n)WT kn .. imWTkn )] n=0 T −1 para k = 0. Esta forma de representación resulta útil. algunas de las propiedades de los números complejos y de las operaciones que pueden ser realizados con los mismos.entonces. T − 1. será necesario realizar operaciones con números complejos. en los próximos acápites marcar se expondrán en forma resumida. . entonces para calcular los diferentes valores. algunos de estos elementos. que implicarán multiplicaciones. se ha empleado la representación de las Series de Fourier. im Xn) (real WTkn .1. . Transformada rápida de Fourier Retomando la expresión que presenta el acápite “Series discretas de Fourier”. pero implica recordar correctamente. A modo de recordatorio. se puede escribir así: síntesis : x ( n) = 1 T T −1 k =0 ∑ x ( k )W − kn T . el par de ecuaciones de la Transformada de Fourier. tendríamos: X (k ) = ∑ [(real ( X (n).. Si descomponemos en partes real (real) e imaginaria (im) la expresión anterior.
Runge (1905) y posteriormente Danielson y Lanczos (1942) describieron algoritmos que lograban que el tiempo de cálculo fuese proporcional a TlogT (Oppenheim A. que ya en 1805. La cantidad de cálculos. que más adelante será analizado. Hemos empleado numerosos de estos algoritmos. C. El período de muestreo de una señal va a determinar un valor que fue objeto de estudio por H.W. aprovechando kn diferentes propiedades de la secuencia WT .W. y por ende del tiempo de cálculo. Un cálculo inicial aproximado. 2 Resolución y “aliasing” La resolución del procedimiento del análisis espectral de una muestra finita (N) de valores de una secuencia temporal dada. entonces el mismo queda definido mediante la expresión: . C. Nyquist (1928) y casi de modo simultáneo por el Académico Kotelnikov en la antigua URSS. y Shafer R. Cooley y J.W.F. Tukey (1965). A partir de los trabajos de J. permite decir que serían necesarias 4T multiplicaciones de números reales y de T(4T – 2) sumas de reales. Gauss se había ocupado en buscar la manera de reducir el número de cálculos a realizar. y actualmente utilizamos un “fabuloso” método descrito por Don Cross (2000). 1989). Si se denomina este valor como Fc. lo que representa un obstáculo serio. que ha sido por ello llamado como frecuencia de Nyquist o frecuencia de Kotelnikov y que constituye la máxima frecuencia que puede ser detectable en el proceso de análisis. se define por la siguiente expresión: Rs = 1 Pm ∗ N donde Rs_resolución y Pm_ período de muestreo. es aproximadamente proporcional al cuadrado de T.V. aún para las computadoras actuales..Desarrollando de modo genérico el producto de números complejos que se muestra. se hizo posible el desarrollo de algoritmos cada vez más eficientes para estos cálculos y que desde entonces han sido denominados como “algoritmos de la transformada rápida de Fourier” (FFT en lengua inglesa). se puede ver que su cálculo sería: [( real X (n) ∗ real WTkn ) − (im X (n) ∗ imWTkn )] + [( real X (n) ∗ imWTkn ) + (im X (n) ∗ real WTkn )] . Está documentado.
el registro original del ECG (continuo) es sometido a discretización para su análisis posterior. distorsionando los resultados. El primero de los casos no es objeto de análisis en esta parte del trabajo.Fc = 1 2 ∗ Pm El teorema del muestreo. o sea. Algo similar ocurre si lo almacenado en la cinta magnética es el valor discretizado de la señal del ECG. discretización de valores de una señal continua en el tiempo.∞ a . entonces la función h(t) estará determinada completamente por sus muestras. si h(f) = 0 para todas las frecuencias mayores que Fc. Si no se toma en cuenta este aspecto del teorema de muestreo. Habitualmente. a partir de la secuencia de valores digitalizados (discretizados). este fenómeno se puede producir si el especialista no presta atención al modo en que fue realizado el registro de la señal original del ECG. A este fenómeno se le denomina en lengua inglesa “aliasing”. se pueden cometer errores importantes por la distorsión que el aliasing produce en los resultados. muestreada con un período Pm. π (t − nPm) Si la función h(t) no está limitada en ancho de banda para los valores menores que la frecuencia de Nyquist (Fc). enunciado por ambos investigadores plantea que una función temporal h(t). pero sí lo será el segundo caso. una serie sometida a discretización. Una secuencia de valores discretos de una serie temporal. El problema del aliasing es una consecuencia del uso de circuitos electrónicos. hay que prestar . tratar de inferir las propiedades de la señal original (continua). Si los sistemas de discretización son aplicados a señales continuas (analógicas). y queremos en la mayoría de los casos. está limitada a un ancho de banda para las frecuencias mayores que Fc. la función h(t) se representa por la expresión: h(t ) = Pm n = −∞ ∑ h( n ) ∞ { sen [2πFc (t − nPm)]} . En este último caso. De forma explícita. por tratarse de la circunstancia con la cual se enfrenta el especialista que estudia la VRC en su trabajo cotidiano. En el caso concreto del análisis de la señal electrocardiográfica (ECG). es de hecho. para el caso de que se trate de señales continuas que muestran el comportamiento de una señal que se produce de manera continua. o sea. y de la introducción del proceso de digitalización de las señales. o cuando se utilizan dispositivos de memoria volátil. la densidad espectral de potencia espectral fuera del rango que va de . las características de la velocidad de la cinta influirán en el proceso de análisis. se desplaza al intervalo entre – Fc y Fc.Fc y de Fc a ∞ . que han sido almacenadas en cinta magnética.
Cuando se emplean los algoritmos de la FFT.81 ms). siempre que ello sea posible. 1990). aunque como antes ya señalamos.. con lo cual se puede comprobar y corregir. si resulta necesario. Existen diversos algoritmos para calcular la duración del ciclo cardíaco. sino que además debe controlarse que el dispositivo de traslación de la cinta. considerado éste como el período de tiempo transcurrido entre el vértice de una onda “R” del complejo “QRS” del electrocardiograma y la inmediata ulterior. el componente N/2 de la serie de valores de salida de la FFT será el que corresponde a la frecuencia de Nyquist. o el reloj interno que regula el proceso de discretización estén funcionando de modo estable. . ya que Fc = 1 2 ∗ Pm y Pm = 1 . el periodo de muestreo debe ser de 2 a 5 milisegundos. puede afectar aún más el efecto del aliasing. tales como: el intervalo P-R. para evitar la distorsión del aliasing en los indicadores espectrales calculados a partir de las mediciones discretas del ECG. Es recomendable. et al. 1990). Veamos algunos ejemplos de interés para el lector. Generalmente...M. Para otras mediciones. Esto podrá ser comprobado en ejemplos que posteriormente serán presentados. cuyos valores son reales. El propio algoritmo utilizado. el Pm de tal proceso debe ser entonces de 0. Fc ∗ 2 En otras palabras. ha sido estudiado para determinar si resultan válidos esos trabajos y sus limitaciones (Merri M.005 Hz.81 ms hayan tenido que aceptarse por razones de orden práctico.01 o 0. En relación con el empleo de diferentes algoritmos para la detección de los vértices de las ondas “R” del ECG. también se han realizado investigaciones (Friesen G. para realizar el análisis de frecuencia de series de cardiointervalos R-R. El hecho de que muchos equipos comerciales utilicen una frecuencia de muestreo de 128 Hz (Pm = 7. cualquier desviación instrumental. De acuerdo con el teorema del muestreo. igualmente resulta importante tener en cuenta el “Pm” empleado para la discretización. que permitirán comprender mejor lo expuesto con anterioridad. como base para el análisis posterior de las series de cardiointervalos R-R.atención no solamente a la frecuencia del proceso de discretización (Pm). emplear como período de muestreo de la señal electrocardiográfica un valor no mayor de 5 ms. se considera que la frecuencia más elevada (rápida) que contiene información útil del ECG es de 50 – 100 Hz. et al.. los estudios realizados con período de muestreo de 7. Un modo comúnmente utilizado es la grabación (en cinta o en memoria) de una señal de tiempo real. intervalos P-P y otros.
(2 ∗ Pm) (2 ∗ 0.25) Puede comprobarse además. Son ellas: Conociendo el límite natural de ancho de banda de la señal a estudiar.25 ∗ 500) 1 1 = = 2 Hz .000976562 0. De acuerdo a lo que ha sido expuesto anteriormente tendremos: Rs = Fc = 1 1 = = 0. Proceder a un filtraje de la señal original. con una frecuencia al . Se desea conocer la resolución que tendría el proceso de análisis espectral de esta serie temporal periódica y la frecuencia de Nyquist correspondiente. se puede decir que existen dos maneras principales de enfrentar la influencia negativa del aliasing. muestrearla con una frecuencia al menos del doble de la frecuencia más alta que la misma posee.08 Hz.008 Hz ∗ 250 = 2 Hz . empleando una frecuencia de discretización de 4 Hz.Si la duración de una serie periódica temporal que se desea someter al análisis espectral tiene una de las duraciones que se muestran en la primera columna de la Tabla 1.0078125 0. Veamos: 0.015625 0.0025 Pasemos ahora a una situación concreta: Se han tomado 500 muestras de una señal sinusoidal de frecuencia 0. ya sea usando métodos de filtraje analógicos o digitales.0050 0. los valores de la resolución del procedimiento serán los que se muestran en la columna correspondiente. Resolución del proceso espectral en dependencia de la duración de la muestra Duración (s) 64 128 256 1024 200 400 Resolución (Hz) 0. Tabla 1.008 Hz . ( Pm ∗ N ) (0. que el componente de frecuencia N/2 (si se aplica la FFT) corresponderá al valor de la Fc. para eliminar las frecuencias indeseadas y proceder luego a realizar un muestreo de la señal filtrada. De manera general.00390625 0.
El tiempo de ejecución de este algoritmo es proporcional a la magnitud del número de estas operaciones que es O(n log(n)). en tanto que una serie de N = 16. que podrían llamarse RealesEnt e ImaginEnt. Ejemplo: el valor real 750 puede ser expresado como el número complejo (750. en una función muestreada en dominio de la frecuencia y que también es compleja. 0). mientras que el número máximo estará en dependencia de las limitaciones de memoria disponible y del sistema operativo del ordenador. El arreglo RealesEnt será el contenedor de los valores reales de la serie temporal que se somete al análisis y el arreglo ImaginEnt. Considerando que en general. Los resultados del procesamiento son entregados al usuario por el algoritmo mediante “n” valores complejos divididos en dos arreglos.000 será inadecuada. Esta limitación está asociada con particularidades del proceso de cálculo que facilita la realización de la FFT y es común a todos sus algoritmos. Algoritmos de la FFT. tiene que ser una potencia entera de 2. los diferentes algoritmos solo varían en cuanto a su implementación y eficiencia. donde O_orden del algoritmo y n_número de elementos en la serie de entrada. El valor de “N”. resulta imprescindible entonces utilizar la propiedad de que un número real puede ser expresado como un número complejo. siendo todos ellos de valor cero (“0”). donde 750 es la parte real y 0 la parte imaginaria de ese número complejo. Algoritmo de Don Cross En un acápite anterior fue señalada la aparición de los algoritmos de la FFT. Como la mayor parte de las veces se opera con valores reales y el algoritmo de la FFT requiere que la señal discretizada esté expresada en valores complejos. deben estar ordenados en dos arreglos4 unidimensionales. que En términos de lenguajes de programación. La FFT es un algoritmo que convierte una función compleja. que representa al número de elementos de la serie compleja temporal. por lo que una serie de N = 13. El algoritmo no permite la entrada de valores que no cumplan esta condición.menos del doble de la frecuencia más alta que pueda quedar después del proceso de filtraje. que tienen que tener el mismo número de elementos (N). 4 .384 será adecuada. vamos a limitarnos en este acápite a la descripción de los rasgos principales del algoritmo desarrollado por Don Cross (2000). El menor número de valores sería “2”. todos del mismo tipo. para su utilización en aplicaciones con el uso de ordenadores. muestreada en tiempo. tendrá el mismo número de elementos. Los valores de la secuencia temporal que se introducen como entrada al algoritmo de la FFT. si como valor imaginario le colocamos el valor cero. los arreglos son estructuras de datos que permiten almacenar una serie de valores. Al finalizar el acápite se incluye el citado algoritmo en su implementación para Pascal.
En la cuarta columna de la tabla se muestra el número de orden de las muestras frecuenciales. aparecen las frecuencias discretas con valores consecutivos múltiplos enteros de la resolución (2Rs.03125 Hz y en las subsiguientes filas de la columna.5 Hz. … . le sigue la que corresponde al valor de la resolución del proceso (“Rs”). con periodo de muestreo de 1 segundo. como son mostrados en la Tabla. que en este caso. lo que resulta del hecho de que son 32 valores discretos (N = 32). así como el cálculo de las densidades espectrales que han sido obtenidos a partir de los datos de salida del algoritmo de la FFT. puede observarse que se han introducido 32 muestras con valores complejos. (2 ∗ Pm) (2 ∗1) En las columnas siguientes a la denominada ImaginEnt. correspondientes a la interpolación de una secuencia de cardiointervalos R-R. como fue antes calculado es de 0. En la Tabla 2 se muestran. Puede observarse que son 32 muestras. para las cuales se han calculado los correspondientes resultados de la FFT. cuyos componentes reales aparecen en la columna con denominación RealesEnt. Puede advertirse que el periodo correspondiente a cada frecuencia . pero la primera frecuencia en aparecer ha recibido el número de orden “0” en lugar de “1”. El número de valores de salida será similar al número de valores de entrada. Por ello. que corresponden a cada valor de frecuencia discreta. el último valor de la columna es el “31” En la quinta columna se muestran los valores de las frecuencias discretas. se han precisado los valores del periodo (1 / Rs). En la sexta columna. a modo de ejemplo concreto. muestreadas con periodo de 1 segundo (“Pm”). los resultados de la aplicación del algoritmo de la FFT (Don Cross) a una serie temporal de 32 valores. mientras que los componentes imaginarios correspondientes. Podemos ver que la primera frecuencia discreta correspondería al valor “0”. los cuales son todos de valor cero (“0”). que corresponden a la salida de la FFT. entonces: Rs = y Fc = 1 1 = = 0. Como la duración total de la serie es de 32 segundos. A los efectos del ejemplo. aparecen mostrados los resultados de la aplicación del algoritmo de la FFT. al igual que las 32 muestras de entrada a la FFT. 3Rs. La secuencia de valores reales de los 32 elementos de la serie temporal se muestra en un histograma secuencial en la Figura 1. aparecen bajo la columna denominada ImaginEnt. nRs).03125 Hz ( Pm ∗ N ) (1∗ 32) 1 1 = = 0.podemos denominar RealesSal e ImaginSal.
210 37.23 1.81981 -25. está teóricamente asociada al valor de la media aritmética de los .671820 -17.68 1.88 1.19549 -1.05764 -81.06250 0.37500 0.42 3.75000 0.1445 -195.42716 -61.0 10.33 1.66 690.39 1. correspondientes a la serie de valores complejos de la serie temporal de entrada.800090 -17.240 563. Resulta evidente que para la primera frecuencia.39399 7.20385 1.96875 1.37 187.28 1.59375 0.48389 135. aparecen los valores de salida del algoritmo de la FFT.50000 0.0000 146.60643 36.19 0. cuyos valores reales e imaginarios se muestran en las columnas segunda y tercera.99338 -105.62500 0.85924 11.52 1.14164 17.4954 -52.45 1.81981 16.17 4.81981 -32.07 1.230 172.00 3.94675 58.87500 0.819810 -11.53107 42.78125 0.94675 36.432830 52.79 1. solamente se obtiene un valor real.90625 0.13 2.00000 -14.54823 -6.42716 -135. el imaginario es cero.28189 -47.81981 53.68 6.56250 0.28125 0.08 2.18750 0.000000 -69.18 2.25000 0.76 323.671820 -10.92056 146.49541 -26.59 22.32 1.58 2.71811 58.57 4. N/O (1) 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 Real Ent (2) 970 980 970 945 965 980 970 930 935 950 960 935 935 955 970 970 945 950 960 960 945 935 935 940 920 915 915 935 920 925 960 980 Imag Ent (3) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 N/O (4) 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 F[i] (5) 0.57 5064.31250 0.discreta se va reduciendo.40625 0.62 13.46875 0.31 48. lo que es una natural consecuencia de la relación inversa entre la frecuencia y el periodo.745570 Imag Sal (8) 0 195.00 1.14 1.10 0. que ocupa el orden “0”.819810 -11.81250 0.14452 11.03 1.46 2.00000 6.56 3.08236 -5.00 Real Sal (7) 30360 42. Esta frecuencia “0” del resultado.28189 -16.10 1.576160 61.71875 0.00000 T (6) 32 16.09375 0.63 4.859240 14.91 2.08236 74.84 12. Tabla 2.68750 0.71811 -105.20385 -17.81981 7.43750 0.53107 53.14164 0.93750 0.570 141.23 1. a medida que se va incrementando el valor correspondiente de la frecuencia.65625 0.800090 -26.67 8.160 20.99338 -52. Ejemplo de aplicación de la FFT a una muestra de 32 valores de una secuencia temporal (Ver explicación detallada en el texto).39399 6.2500 Amplitud Espectro (10) 948.28 2.44 En las columnas séptima y octava.00 5.60643 250.21875 0.0000 -81.18 1.53125 0.84375 0.64 0.74557 250.19549 6.54823 5.33 4.81981 49.60 1.34375 0.05764 -25.64 2.24 1643.460 45.43283 69.03125 0.92056 49.57616 -74.4834 Potencia Espectro (9) 57608100 2502.000000 47.19 7.20 2.75 8.12500 0.67 2.78 1.
Debe también ser objeto de atención. tanto reales como imaginarios. Más adelante volveremos sobre esta cuestión. fue expuesto que dos números complejos son conjugados si se da la condición ( a .valores de la serie de entrada. Su magnitud (amplitud) es. o valor de la fase correspondiente). correspondiente a la frecuencia “N/2”. constituyendo una imagen en espejo de los correspondientes a la primera mitad. podemos reducir su valor para una serie temporal de entrada dada. resulta opuesto al que tuvo en la primera mitad. que el valor del signo correspondiente a los valores imaginarios en la segunda mitad de la serie. resulta entonces que las salidas para las frecuencias positivas y negativas son similares y podemos decir que redundantes. respecto a los signos en la segunda mitad. si a dicha serie se le calcula la media aritmética y luego se resta la misma de cada valor de la serie original. − b ) . puede irse resumiendo que los valores a la salida del algoritmo de la FFT muestran una imagen especular de valores. cuya longitud sea diferente de 2N. Resulta. Vale recordar que en la parte que fue dedicada a la actualización acerca de la teoría de los números complejos. . También es posible. b ) y ( a. lo que significa que sus partes reales son similares y que sus partes imaginarias tienen signo contrario. cuando los valores de entrada han sido reales y que solo se diferencian en el signo. Por tanto. usando un procedimiento especial del algoritmo. que estos valores son números conjugados complejos. Ello es consecuencia de que los valores de entrada a la FFT. que recordemos algunas expresiones de los números complejos que nos serán de utilidad inmediata. Una inspección de los valores calculados. fueron reales. en todo caso. se repiten en la segunda mitad de la serie de valores en la salida de la FFT. Como en la inmensa mayoría de los casos emplearemos valores reales a la entrada de la FFT. calcular la FFT de secuencias de números reales enteros. Nos parece conveniente ahora. resultando que los signos de los valores imaginarios son opuestos en la segunda mitad. b) = ρ ϕ . muestra que las cifras obtenidas a partir de la frecuencia discreta “N/2”. El algoritmo de Don Cross. Podemos considerar el valor real e imaginario de salida de la FFT para cada frecuencia discreta en la forma “modulo-argumental” que se expuso en acápites anteriores: ( a. es cero. proporcional a la media aritmética de los valores de entrada como ya se dijo. el hecho de que el valor de salida imaginario. y por ello. También es de notar. tiene algunas facilidades que permiten que los cálculos se limiten a un solo componente de frecuencia discreta que se desee evaluar (ya sean valores de amplitud de componentes real e imaginario.
a_ parte real. bastaría con calcularlo a partir de cualquiera de las expresiones: . calculando la raíz cuadrada se obtendría ρ.3389963 = ρ = ± 2 775. elevado al cuadrado. Consideremos como ejemplo los valores real e imaginario correspondientes a la salida de la FFT para la frecuencia discreta 0.Figura 1 Representación mediante un histograma secuencial de los valores de la secuencia utilizada en el texto como un ejemplo. correspondería al valor imaginario.49541) 2 + (11. sumado con el cuadrado del valor imaginario y después de ser sumados ambos. Para encontrar el argumento ϕ.3750 en la Tabla 2.19549) 2 = ρ = ± 2 650. Para obtener el módulo del resultado basta con operar: ρ = ± 2 (−25.0159311 + 125. se desprende que el módulo ρ estará dado por el valor real obtenido para la frecuencia discreta correspondiente. En este caso. correspondería al valor real obtenido y b_ parte imaginaria.3549274 = ± 27.84519577 . Recordando que ρ = 2 a2 + b2 .
En el ejemplo que venimos analizando en este acápite. El punto ρϕ en un plano de coordenadas. en tanto que las de frecuencia-fase-intensidad. Los valores principales de la función “arc tan” tienen que cumplir la condición: − π 2 p arc tan x p π 2 . resultan absolutamente equivalentes a la representación en coordenadas de frecuencia-fase-intensidad. cuál es la información de fase. el afijo se encontrará en el segundo cuadrante. posibilidad que también brinda el algoritmo de Don Cross que venimos analizando. En el caso de la frecuencia discreta.= − 0.19549. Debe advertirse. aplicando el proceso de la FFT inversa.31 . para la cual se efectuaron algunos cálculos anteriormente. todos aquellos algoritmos que usan el concepto de la FFT. no queda claro de la simple inspección de la tabla. Como hemos visto en otro acápite anterior. fue mencionado el caso de la transformada de Fourier. este valor ρϕ son las coordenadas polares del número complejo correspondiente. Se indicó que los valores de la serie temporal. es el afijo que corresponde en este caso a la frecuencia discreta en cuestión. que a los efectos del cálculo de la fase. ese valor principal de la fase correspondería a – 0.43912 : a ρ sen ϕ = b . Los valores reales e imaginarios del valor complejo de salida de la FFT para cada frecuencia discreta son portadores.49541 y en el de las ordenadas de 11. representados en las coordenadas tiempo-intensidad. . En este caso. sin embargo. con valor en el eje de las abscisas de – 25. ρ cos ϕ = a . Sin embargo. lo están en las columnas séptima y octava. expresado generalmente en radianes y obtenido mediante la expresión ya antes mostrada: tan ϕ = b ϕ = arctan b a . y en general. Teniendo estos elementos. ϕ = arc tan b a = 156 o 17 ' 34 ". por lo tanto. al señalar un ejemplo del primer tipo de los mismos. La fase asociada a cada valor de frecuencia discreta está dada por el valor del ángulo ϕ.414 radianes en la frecuencia discreta 11. En el acápite “Métodos de Análisis de Procesos”. se puede reconstituir unívocamente la serie original de valores de entrada a la FFT. se toman solo los denominados valores principales de la función “arco tangente”. de la información referente a la intensidad y a la fase. las coordenadas tiempo-intensidad de la serie de entrada están representados por los que se muestran en las columnas segunda y tercera de la Tabla 2.
048 0.355 0.175 0.578 -0. Por ello.122 0. En las abscisas se han situado los valores de las frecuencias discretas (Hz). se muestra una curva que describe el perfil de la variación de la fase de la secuencia de entrada. Este diagrama se denomina diagrama de fase y .440 0. se obtienen los valores que se muestran tabulados a la derecha de la Figura 2.414 0.068 1. En la parte izquierda.901 0.210 -1.131 -0.495 0. calculando la fase para las 15 primeras frecuencias discretas del ejemplo.1.782 -1. en tanto en las ordenadas. aparecen los valores de la fase expresados en radianes.033 -0.430 Figura 2 Diagrama de fases del espectro de la secuencia utilizada en el ejemplo. Figura 3 Diagrama de afijos de la serie de salida de la FFT del ejemplo en el texto.
var ImagOut: array of double ). var RealOut: array of double.we violate array bounds rules *) unit Fourier. interface (*-----------------------------------------------------------------procedure fft Calculates the Fast Fourier Transform of the array of complex numbers represented by 'RealIn' and 'ImagIn' to produce the output complex numbers in 'RealOut' and 'ImagOut'.junto al espectro de potencia.. ------------------------------------------------------------------*) procedure fft_integer ( NumSamples: word. En la Figura 3 se muestra un diagrama de afijos que se corresponde con los valores calculados en el ejemplo. (*-----------------------------------------------------------------procedure fft_integer Same as procedure fft. var ImagIn: array of integer. but uses integer input arrays instead of double..Don Cross <dcross@intersrv. var ImagOut: array of double ). Make sure you call fft_integer_cleanup after the last time you call fft_integer to free up memory it allocates. { must be a positive integer power of 2 } var RealIn: array of double. veremos la forma de calcular el espectro de potencia. var ImagOut: array of double ). . var ImagIn: array of double. var RealOut: array of double. var RealIn: array of integer.com> ==================================================================*) {$N+. Más adelante.E+} (* Allows code to use type 'double' and run on any iX86 machine *) {$R-}(* Turn off range checking. (*-----------------------------------------------------------------procedure ifft Calculates the Inverse Fast Fourier Transform of the array of complex numbers represented by 'RealIn' and 'ImagIn' to produce the output complex numbers in 'RealOut' and 'ImagOut'. { must be a positive integer power of 2 } var RealIn: array of double.pas . var RealOut: array of double. constituyen los resultados principales del análisis espectral. -----------------------------------------------------------------*) procedure fft ( NumSamples: word. var ImagIn: array of double. Algoritmo de Don Cross (Código) fourier. -----------------------------------------------------------------*) procedure ifft ( NumSamples: word.
Use this instead of 'fft' when you only need one or two frequency samples. var i. you could calculate the DFT of 100 points instead of rounding up to 128 and padding the extra 28 array slots with zeroes. y: word. (*-----------------------------------------------------------------procedure CalcFrequency This procedure calculates the complex frequency sample at a given index directly. var ImagOut: double ). end. not the whole spectrum. ------------------------------------------------------------------*) procedure fft_integer_cleanup. { can be any positive integer } FrequencyIndex: word. end. begin y := 2. . ------------------------------------------------------------------*) procedure CalcFrequency ( NumSamples: word. end. y := y SHL 1. ReverseBits := rev. rev: word. function ReverseBits ( index. begin for i := 0 to 16 do begin if (PowerOfTwo AND (1 SHL i)) <> 0 then begin NumberOfBitsNeeded := i. var ImagIn: array of double. you must call 'fft_integer_cleanup' after the last time you call 'fft_integer' in order to free up dynamic memory. end. exit. var i: word. For example. exit. var RealOut: double. end. NumBits: word ): word. end. var i. for i := 1 to 15 do begin if x = y then begin IsPowerOfTwo := TRUE. It is also useful for calculating the Discrete Fourier Transform (DFT) of a number of data which is not an integer power of 2.. begin rev := 0. end. NumSamples-1} var RealIn: array of double. implementation function IsPowerOfTwo ( x: word ): boolean. for i := 0 to NumBits-1 do begin rev := (rev SHL 1) OR (index AND 1). function NumberOfBitsNeeded ( PowerOfTwo: word ): word. { must be in the range 0 .(*-----------------------------------------------------------------procedure fft_integer_cleanup If you call the procedure 'fft_integer'. IsPowerOfTwo := FALSE. index := index SHR 1.
(alpha*ai . i := 0. NumBits := NumberOfBitsNeeded (NumSamples). procedure FourierTransform ( AngleNumerator: double. var ImagIn: array of double. var RealOut: array of double. while BlockSize <= NumSamples do begin delta_angle := AngleNumerator / BlockSize.0 * alpha * alpha. RealOut[j] := RealIn[i]. alpha := sin ( 0. ai := ai . ar := ar . ImagOut[j] := ImagOut[j] + ti. BlockSize. BlockEnd := BlockSize. while i < NumSamples do begin ar := 1. end. delta_ar := alpha*ar + beta*ai. NumBits ). BlockSize := BlockSize SHL 1. delta_ar: double. ti := ar*ImagOut[k] + ai*RealOut[k].' ). for i := 0 to NumSamples-1 do begin j := ReverseBits ( i. end. beta: double. halt. var NumBits. begin if not IsPowerOfTwo(NumSamples) or (NumSamples<2) then begin write ( 'Error in procedure Fourier: NumSamples='. ai: double. NumSamples: word.0. end. beta := sin ( delta_angle ). ImagOut[k] := ImagOut[j] . writeln ( ' is not a positive integer power of 2. alpha. i := i + BlockSize. .end. (* sin(0) *) j := i. RealOut[j] := RealOut[j] + tr. (* cos(0) *) ai := 0. tr.0.ti. end.beta*ar). delta_angle. ar. RealOut[k] := RealOut[j] . INC(j).5 * delta_angle ). end.tr.delta_ar. BlockEnd := 1. var RealIn: array of double.ai*ImagOut[k]. procedure fft ( NumSamples: word. BlockSize := 2. j. ti. var RealIn: array of double. var ImagIn: array of double. alpha := 2. end. ImagOut[j] := ImagIn[i]. var ImagOut: array of double ). n. NumSamples ). i. for n := 0 to BlockEnd-1 do begin k := j + BlockEnd. tr := ar*RealOut[k] . k. BlockEnd: word.
TempArraySize * sizeof(double) ). RealOut. ImagTemp^[i] := ImagIn[i]. var i: word. var ImagIn: array of integer. ImagOut ). procedure ifft ( NumSamples: word. NumSamples. var i: word. for i := 0 to NumSamples-1 do begin RealTemp^[i] := RealIn[i]. RealOut. end. begin if TempArraySize > 0 then begin if RealTemp <> NIL then begin FreeMem ( RealTemp. end. GetMem ( ImagTemp.. var RealIn: array of integer. var RealOut: array of double. var RealIn: array of double. TempArraySize := NumSamples. procedure fft_integer_cleanup. NumSamples. ImagOut ). TempArraySize: word. } GetMem ( RealTemp. RealTemp^. { free up memory in case we already have some. var RealOut: array of double. ImagIn. NumSamples * sizeof(double) ).0] of double. var RealTemp. .var RealOut: array of double.. end. NumSamples * sizeof(double) ). var ImagOut: array of double ). RealOut. ImagOut[i] := ImagOut[i] / NumSamples. begin FourierTransform ( 2*PI. var ImagOut: array of double ). ImagTemp^. RealIn. ImagTemp: ^doubleArray. end. type doubleArray = array [0. FourierTransform ( 2*PI.. ImagIn. (* Normalize the resulting time samples. procedure fft_integer ( NumSamples: word. end. end. *) for i := 0 to NumSamples-1 do begin RealOut[i] := RealOut[i] / NumSamples. var ImagIn: array of double. begin if NumSamples > TempArraySize then begin fft_integer_cleanup. RealIn. ImagOut ). begin FourierTransform ( -2*PI. NumSamples. RealTemp := NIL. var ImagOut: array of double ).
for k := 0 to NumSamples-1 do begin { Update trig values } sin3 := beta*sin2 . sin2 := sin ( -theta ). cos3 := beta*cos2 . end. end. procedure CalcFrequency ( NumSamples: word. beta := 2 * cos2. {flag that buffers RealTemp. end. NumSamples-1 } var RealIn: array of double.ImagIn[k]*sin3. sin1 := sin2. begin { Unit initialization code } TempArraySize := 0. var ImagIn: array of double.cos1. var k: word. var ImagOut: double ). cos2 := cos3. cos3. cos2 := cos ( -theta ). TempArraySize := 0. ImagTemp := NIL. end. beta: double. cos2. sin2 := sin3. TempArraySize * sizeof(double) ). RealOut := RealOut + RealIn[k]*cos3 . cos1 := cos2. theta := 2*PI * FrequencyIndex / NumSamples. sin2. sin3: double. ImagTemp := NIL. { must be in the range 0 . cos1 := cos ( -2 * theta ). RealImag not allocated} RealTemp := NIL.end. ImagOut := 0.0. if ImagTemp <> NIL then begin FreeMem ( ImagTemp. ImagOut := ImagOut + ImagIn[k]*cos3 + RealIn[k]*sin3. end. { must be integer power of 2 } FrequencyIndex: word. end. begin RealOut := 0. var RealOut: double. sin1 := sin ( -2 * theta ). cos1..0. sin1. .sin1. theta.
Estévez Báez M.. Shafer R. IEEE Trans Biomed Eng 1990. 5. 6.. Denisova V. Movsisyants S. 4. 37:99–106. and Tukey J. Merri M.I. (1989) Discrete-time signal processing. Titlebaum EL.:Prentice-Hall. Friesen GM. Quint SR. 5.8. 2. and Tukey J. 37:85– 98. URSS. (1965) An algorhythm for the machine computation of complex Fourier series.W. Yates SL.B.W. Jannett TC. 3. 851-860 (en idioma ruso). New York. Nikitina L.W. . (1982) Particularidades de los procesos extralentos en algunas enfermedades del cerebro.V.A.. N. Oppenheim A.J. Sampling frequency of the electrocardiogram for the spectral analysis of heart rate variability. Farden DC. Dover Publications.. (1958) The measurement of power spectra. Englewood Cliffs. A comparison of the noise sensitivity of nine QRS detection algorithms. Cooley J. Mathematics of Computation Vol 19:297-301. Blackman R. No. Nogle HT. IEEE Trans Biomed Eng 1990.V.Bibliografía 1. pp. Mottley JG. Jadalloh MA. Revista de Fisiología Humana. T.W.
Documents Similar To AnalisisEspectralMedianteUsoFFT
sofia aillapan
More From lucho999

References: Resolución 
 resolución 
 resolución 
 resolución 
 Resolución 
 Resolución 
 resolución 
 resolución