AMPLIACIÓN DE MATEMÁTICAS (2o Ing. de Telecomunicación) Departamento de Matemática Aplicada II. Universidad de Sevilla CURSO ACADÉMICO 2008-2009

Práctica IV: La transformada Discreta/Rápida de Fourier.

Interpolación Trigonométrica Compleja Interpolación trigonométrica compleja. El problema de la interpolación trigonométrica compleja, con periodo prefijado T > 0, consiste en, dado un vector arbitrario y = (y0, y1 , . . . , yN −1 )t ∈ CN , encontrar una combinación lineal de exponenciales complejas ¶ ¶¶ µ µ µ 2π 2π β0 + β1 exp ti + . . . + βN−1 exp (N − 1)ti T T ¶ µ N −1 2π 1X βn exp n ti = N n=0 T

1 P (t) = N

tal que P

µ

k T N

¶

= yk ,

para todo k = 0, 1, . . . , N − 1.

Conviene mencionar que el factor inicial 1/N no es más que un convenio y se incluye para seguir la notación habitual en los textos de procesado digital de señales; convenio que, por otra parte, también utiliza el programa MATLAB. Es bien conocido que este problema de interpolación tiene siempre solución única y ello se debe a las notables propiedades de las matrices de Fourier. ⎤ ⎡ 1 1 1 ··· 1 ⎥ ⎢ 1 ω ω2 ... ωN −1 ⎥ ⎢ 4 2(N−1) ⎥ ⎢ 1 ω2 ω · · · ω FN := ⎢ ⎥. ⎥ ⎢ .. .. .. .. . . ⎦ ⎣ . . . . . N−1 2(N−1) (N−1)(N −1) ω ··· ω 1 ω 1

¶ 2πi , es decir, la raíz primitiva N-ésima de la unidad. En Como es usual, ω = exp N concreto, se verifica que FN−1 = NFN e µ

y=

1 FN β, N

β = FN y,

donde β = (β0 , β1 , . . . , βN−1 )t . Ejercicio resuelto 1. Verifique numéricamente que F4 F4 = 4id4 . >> fou4=[1 1 1 1;1 i -1 -i;1 -1 1 -1;1 -i -1 i] >> fou4*conj(fou4) Ejercicio 2. Repita el ejercicio anterior con F6 . La transformada rápida de Fourier. Acabamos de ver que la resolución efectiva del problema de la interpolación trigonométrica compleja puede realizarse simplemente multiplicando por la matriz de Fourier conjugada, lo cual conduce a N 2 mulpiplicaciones complejas y otras tantas sumas. Para las aplicaciones actuales este orden computacional es realmente alto pero, afortunadamente, hay un algoritmo que es capaz de reducir este costo a un nivel de aproximadamente N log2 N: la transformada rápida de Fourier. La orden fft de MATLAB implementa este método numérico. Ejercicio resuelto 3. Resuelva el problema de interpolación trigonométrica compleja en [0, 2π] para el vector y = (1, 2, 3, 4)t y compruebe gráficamente que la parte real del correspondiente "polinomio" de interpolación complejo interpola los pares ½ ¾ π 3π (0, 1), ( , 2), (π, 3), ( , 4) . 2 2 ¿Qué pares interpola la parte imaginaria?

(En un fichero de datos intercom.m) y=1:4;beta=fft(y); t=0:0.001:2*pi; fun=zeros(1,length(t)); for k=1:4 fun=fun+beta(k)*exp((k-1)*t*i); end fun=0.25*fun; plot(t,real(fun)) hold on plot(0:pi/2:3*pi/2,1:4,’or’) hold off shg pause plot(t,imag(fun)) 2

hold on plot(0:pi/2:3*pi/2,zeros(1,4),’or’) hold off shg

Interpolación Trigonométrica Real Interpolación trigonométrica real. Sea f : R −→ R una señal periódica de periodo T > 0, de la cual sólo se conocen sus valores en una cierta colección equiespaciada de puntos k T, k = 0, 1, . . . , N − 1. N Por comodidad, asumiremos siempre que N = 2M es par. Por otro lado, consideremos el conjunto Trig(T, N) de las siguientes funciones T -periódicas tk :=

½ µ ¶ µ ¶¾ ½ µ ¶ µ ¶ ¾ 2π 2π 2π 2π cos 0 t = 1, cos M t ∪ sen n t , cos n t : n = 1, . . . , M − 1 . T T T T El problema de la interpolación trigonométrica real consiste en encontrar una combinación lineal de las funciones del conjunto Trig(T, N) (denominada polinomio trigonométrico equilibrado)

P (t) = a0 +

M−1 Xµ n=1

¶ ¶¶ ¶ µ µ µ 2π 2π 2π an cos n t + bn sen n t + aM cos M t , T T T

tal que P (tk ) = f (tk ),

para todo k = 0, 1, . . . , N − 1.

Las condiciones de interpolación anteriores, además, muestran claramente el por¡ qué ¢de 2π no incluir ¡ 2π en ¢ el conjunto Trig(T, N) las funciones elementales y T -periódicas sen 0 T t y sen M T t . Simplemente son funciones redundantes en el proceso de interpolación pues valen 0 en cada tk . Desde el punto de vista del Análisis de Fourier, lo que se pretende con la interpolación trigonométrica es conseguir una “serie finita trigonométrica” fácilmente manipulable, que represente a f y que pueda verse, por tanto, como una versión truncada de la serie de Fourier de f en [0, T ]. En cualquier caso, conviene subrayar que, en general, esta serie “finita trigonométrica” no es necesariamente una suma parcial finita de la serie de Fourier asociada a f . Puede probarse que el polinomio trigonométrico equilibrado P siempre existe y es único, o equivalentemente, los N coeficientes reales

3

Γ := {ao , a1 , . . . , aM , b1 , . . . , bM−1 } están unívocamente determinados. Fórmulas de interpolación usando la orden fft. Por comodidad, supondremos otra vez que T = 2π y consideremos la correspondiente función f : R → R periódica de periodo 2π. A su vez, denotemos µ ¶ k yk := f 2π , k = 0, 1, . . . , N − 1 (N = 2M). N

Si escribimos y = (y0, y1 , . . . , yN−1 )t ∈ RN , podemos considerar su transformada discreta de Fourier β = f f t(y) ∈ CN y obtener, por tanto, el polinomio de interpolación trigonométrico complejo asociado al vector y. En concreto, 1X h(t) := βn exp(nti). N n=0 N−1

Al ser y un vector real, puede probarse que las coordenadas de β verifican que β0 , βM ∈ R y βn = βN−n , para n = 1, . . . , M − 1. La utilidad de este aparente rodeo es que a partir de β es inmediato el cálculo de los coeficientes del problema de interpolación trigonométrica real. A saber

a0 :=

β0 , N

an :=

2 Re(βn ) , N

bn :=

−2 Im(βn ) N

(n = 1, . . . , M − 1),

aM :=

βM , N

donde, lógicamente, {a0 , a1 , . . . , aM , b1 , . . . , bM−1 } son los coeficientes del polinomio de interpolación trigonométrico real asociado a y en [0, 2π]. Ejercicio resuelto 4. Diseñe una función en MATLAB que implemente la interpolación trigonométrica real a partir de la transformada rápida de Fourier. La función debe determinar, al menos, los coeficientes del polinomio de interpolación trigonométrico equilibrado y una gráfica de dicho polinomio en [0, T ], donde T es el periodo de muestreo. Los argumentos de entrada deben ser la muestra, el periodo de muestreo T y el paso h del mallado para generar la gráfica del polinomio en [0, T ]. Usando la orden plot, represente gráficamente, en el intervalo [0, 0.32] el siguiente polinomio trigonométrico de periodo 0.32 = 1/f1 : f(t) = 0.5 + 2 sen(2πf1 t) + cos(2πf2 t),

f1 = 3.125Hz, f2 = 6.25Hz.

A partir de f , obtenga una muestra de 64 puntos equiespaciados en el intervalo anterior. Verifique que la función diseñada permite reconstruir gráficamente, a partir de la muestra, la función f. (Gráfica de la señal original) >> t=0:0.001:0.32; 4

>> f1=3.125;f2=6.25; >> y=0.5 +2*sin(2*pi*f1*t)+cos(2*pi*f2*t); >> plot(t,y) (Obtención de la muestra) >> t=0:0.32/64:0.32;t=t(1:64); >> m=0.5+2*sin(2*pi*f1*t)+cos(2*pi*f2*t); (Fichero para la interpolación trigonométrica real) function [a,b]=coeftrig(x,T,h) N=length(x);M=N/2; beta=fft(x); a(1)=beta(1)/N; a(M+1)=beta(M+1)/N; for i=2:M a(i)=2*real(beta(i))/N; b(i)=-2*imag(beta(i))/N; end %%%%%%%%%%%Dibujo del polinomio%%%%%%%%%%% t=0:h:T; p=2*pi/T; L=T/h +1; u=zeros(1,L); for k=1:M-1 u=u+a(k+1)*cos(p*k*t)+b(k+1)*sin(p*k*t); end u=u+a(1)+a(M+1)*cos(p*M*t); b=b(2:M); plot(t,u),shg Nota: Utilice la función con los siguientes valores: x = m, T = 0.32 y h = 0.001. Ejercicio resuelto 5. Considere nuevamente la señal f(t) = 0.5 + 2 sen(2πf1 t) + cos(2πf2 t),

f1 = 3.125Hz, f2 = 6.25Hz.

Con la muestra de 64 puntos equiespaciados obtenida en el ejercicio anterior, realice un periodograma (gráfica de amplitud frente a frecuencia) del correspondiente polinomio de interpolación trigonométrico. Compruebe la exactitud del análisis espectral, es decir, que en la señal reconstruida aparecen las mismas frecuencias que en la original. function barfou(x,v) N=length(x);M=N/2; beta=fft(x); beta(1)=beta(1)/2; beta(M+1)=beta(M+1)/2; beta=(2/N)*beta; 5

beta=beta(v); bar(real(beta),’b’); hold on bar(imag(-beta),’r’); hold off shg Nota: Ejecute el fichero anterior con v=1:32 y v=1:64. Ejercicio 6. Realice un análisis espectral de una onda triangular de amplitud ±1 y periodo un segundo, muestreada a 1/32 en un ciclo. Compare los resultados obtenidos en el desarrollo en serie de Fourier de dicha señal. ¿Hay discrepancias significativas?

Análisis Espectral de Señales Análisis espectral de señales. La transformada discreta/rápida de Fourier se usa de modo masivo en multitud de temas relacionados con el procesamiento digital de señales analógicas. De obligada mención son los análisis y síntesis espectrales de señales, la correlación cruzada de señales o la convolución de señales. Recordemos que una señal analógica es una función continua del tiempo t ∈ R −→f(t) ∈ R que representa información, como el sonido de una voz, la presión sanguinea,... Para procesar esta información con un computador, se toma una muestra de la señal cada T segundos y así se genera una cierta señal digitalizada. Puesto que tomamos muestras cada T segundos, hay T −1 muestras por segundo y, se dice, entonces que la frecuencia de muestreo es de T −1 Hz. En la práctica, puede asumirse que las señales más utilizadas son las aperiódicas de banda limitada y las periódicas finitas, es decir, formadas por un número finito de armónicos. Si la correspondiente muestra consta de un total de N valores, entonces el n-ésimo valor es yn = f(nT ),

n = 0, 1, · · · , N − 1.

La transformada discreta/rápida de Fourier nos permite convertir la señal digital anterior (yn ) en el dominio del tiempo en un conjunto de puntos (βn ) que representan el contenido en frecuencia. Puesto que los puntos en los que tomamos las muestras están igualmente espaciados en el intervalo temporal [0, NT ], los coeficientes calculados con la transformada discreta correspondiente a frecuencias separadas por (NT )−1 Hz. Fenómenos específicos: aliasing, leakage. La frecuencia de Nyquist de una señal es el doble de la máxima frecuencia presente en dicha señal. El teorema del muestreo nos dice que para muestrear “adecuadamente” una señal debemos muestrear al menos con la frecuencia de Nyquist. Cuando se muestrea por debajo de la frecuencia de Nyquist se produce el fenómeno denominado aliasing (palabra que procede de “alias”). Este fenómeno consiste en que las frecuencias superiores a la mitad de la frecuencia de muestreo aparecen como frecuencias distintas en la senal muestreada (esto es con “un alias”). 6

Supongamos que tenemos una sunisoide de 10Hz con una duración de 10 segundos y que la muestreamos de manera muy fina, por ejemplo, 1000Hz : >> sf=1000; t=0:1/sf:1; >> f=10; y=sin(2*pi*f*t); Si diezmamos la senal discreta “y” tomando una muestra de cada k tenemos la secuencia >> z=y(1:k:end); La frecuencia de muestreo de la secuencia z será k veces menor que la de y : fz =

1000 fy = Hz. k k

Puesto que la señal original tiene 10Hz de frecuencia, la menor frecuencia de muestreo aconsejable será de 20Hz, es decir, que como máximo podemos tomar k = 50. Si muestreamos tomando k > 50 aparecerá aliasing. Otro aspecto a tener en cuenta al muestrear es el denominado fenómeno de leakage. Si f es periódica de periodo p, el problema surge cuando T no es múltiplo de p. En concreto, cuando no se muestrea sobre un número entero de ciclos de los armónicos presentes en la señal, aquellas componentes de la DFT más cercanas a dichos armónicos se “dispersan o derraman” a otras frecuencias. La razón del fenómeno estriba en que, al forzar en el muestreo un periodo distinto del propio de la señal, intentamos realmente reconstruir una señal más compleja (usualmente con más discontinuidades) que la original y que requiere, por tanto, de un mayor número y variedad de frecuencias. El leakage en el caso aperiódico surge, simplemente, de que tenemos necesariamente que muestrear en un intervalo finito [0, L = NT ] y la señal puede tomar valores no nulos en toda la recta real o en un intervalo muy grande. Nuevamente, esto puede forzar discontinuidades no presentes en la señal y, de ahí, nuevamente el “derrame” a otras frecuencias. Un modo de paliar este fenómeno es utilizar funciones (window functions) ω(t) en [0, L] que decaigan suavemente en los extremos del intervalo y considerar, como muestras de la señal, los valores (ω(nT )f(nT )). Ejercicio resuelto 7. Determine la transformada discreta/rápida de Fourier de 512 puntos muestreados sobre un periodo de un segundo y procedentes de la señal f (t) = sen(2πf1 t) + 2 sen(2πf2 t), donde f1 = 30Hz y f2 = 400Hz. Describa el fenómeno de aliasing que aparece y como lo podría solventar. (Gráfica de la señal original) >> t=0:0.001:1; >> f1=30;f2=400; >> y=sin(2*pi*f1*t)+2*cos(2*pi*f2*t); >> plot(t,y) Nota: La frecuencia de 400Hz hace que la visualización en el tiempo sea confusa. (Obtención de la muestra: 512 puntos) >> t=0:1/512:1;t=t(1:512); 7

>> m=sin(2*pi*f1*t)+2*sin(2*pi*f2*t); (Análisis espectral) barfou(m,1:206) Nota: A tenor del análisis espectral, la “señal reconstruida” es g(t) = sen(2πf1 t) − 2sen(2πg2 t),

con g2 = 212.

Hay un fenómeno de aliasing debido a un muestreo lento. Para solventar el problema basta muestrear con 802 puntos puesto que la frecuencia más alta es de 400Hz. Ejercicio 8. Considere ahora una muestra de 256 puntos equiespaciados de la señal h(z) = sen(2πf3 t),

f3 = 30.27Hz

en el intervalo [0, 1] y obtenga el correspondiente periodograma. Determine cúal o cuales de las componentes mostradas deben considerarse erróneas y explique razonadamente por qué aparecen.

8