Análisis de Fourier

Vamos a ver cómo

  • Serie de Fourier de una señal acotada (o periódica)
  • Transformada de Fourier
  • Fast Fourier Transform
  • Teorema de Nyquist, reconstrucción a partir de una muestra equiespaciada
  • Aplicación a oleaje aleatorio

Los métodos de cuadratura que hemos visto aproximan la integral de una función de variable real, y a valores reales. Al usar la transformada de Fourier, aparecen funciones de variable real, pero a valores complejos. Integramos por separado la parte real y la parte imaginaria...

In [22]:
function integral = complex_quadrature(func, a, b)
    real_func = @(x)(real(func(x)));
    imag_func = @(x)(imag(func(x)));
    [real_integral, error_real] = quad(real_func, a, b);
    [imag_integral, error_imag] = quad(imag_func, a, b);
    integral = real_integral + 1j*imag_integral;
end

Serie de Fourier de una función acotada, o periódica

Comenzamos por una integral de una función en un intervalo [-L,L]. Podemos asumir que la función se extiende a todo $\mathbb{R}$ de forma periódica, y construímos algunos términos de su serie de Fourier: $$ \hat{x}_k = \frac{1}{L}\int_{0}^L x(t)e^{-i2\pi \frac{k}{L} t}\:dt $$ La integral se puede hacer entre -L/2 y L/2 porque la función $x$ es periódica. En condiciones bastante generales podemos reconstruir la función $x$ de forma perfecta: $$ x(t) = \sum_{k=-\infty}^\infty \hat{x}_k e^{i2\pi \frac{k}{L} t} $$ y si sumamos unos cuantos términos, reconstruimos $x$ de forma aproximada: $$ x(t) = \sum_{k=-N}^N \hat{x}_k e^{i2\pi \frac{k}{L} t} $$

Ejemplo

In [23]:
L = 2*pi;
f = @(x)(1 + exp(-x.^2).*sin(5.*x));

t_eval = linspace(-L/2, L/2, 200);
f_eval = f(t_eval);
plot(t_eval,f_eval)

Para empezar, vamos a calcular las integrales que definen los coeficientes $\hat{x}_k$ usando los métodos de cuadratura quad recomendados por nuestras librerías.

Ejercicio

  • Comprueba cómo mejora la aproximación al aumentar el número de términos en la aproximación mediante la serie de Fourier.
In [25]:
orden = 6;
ks = [-orden:orden];
Xk = zeros(2*orden+1);
integrando = f;
integral = complex_quadrature(integrando,-L/2,L/2);
Xk(1) = (1/L)*integral;
for k = ks
    integrando = @(t) (f(t).*exp(-1j*(2*pi*k/L)*t));
    integral = complex_quadrature(integrando,-L/2,L/2);
    if k<0
        Xk(2*orden+2+k) = (1/L)*integral;
    else
        Xk(1+k) = (1/L)*integral;
    end
end

hold on;
Ntimes = 200;
t_eval = linspace(-L/2, L/2, Ntimes);
f_eval = f(t_eval);
plot(t_eval,f_eval)

%reconstruccion
fourier = zeros(1, Ntimes);
for k=ks
    onda = @(t)(exp(1j*(2*pi*k/L)*t));
    onda_eval = onda(t_eval);
    if k<0
        fourier = fourier + Xk(2*orden+2+k)*onda_eval;
    else
        fourier = fourier + Xk(1+k)*onda_eval;
    end
end
legends = cell(1, 2);
legends{1} = 'f';
legends{2} = sprintf('Serie de Fourier hasta grado %d',orden);
max(abs(imag(fourier)))%comprobamos que la parte imaginaria es cero
plot(t_eval, real(fourier), 'g-')
title('Reconstruccion con la serie de Fourier')
legend(legends)
hold off;

ans =    4.4235e-17

Fenómeno de Gibbs

Si la función es discontinua, aparece el fenómeno de Gibbs.

Ejercicio

  • Comprueba el fenómeno de Gibbs al aumentar el número de términos en la aproximación mediante le serie de Fourier.
In [4]:
L = 2*pi;
f = @(x)(1 + sign(sin(2*x)));

orden = 20;

ks = [-orden:orden];
Xk = zeros(2*orden+1);
integrando = f;
integral = complex_quadrature(integrando,-L/2,L/2);
Xk(1) = (1/L)*integral;
for k = ks
    integrando = @(t) (f(t).*exp(-1j*(2*pi*k/L)*t));
    integral = complex_quadrature(integrando,-L/2,L/2);
    if k<0
        Xk(2*orden+2+k) = (1/L)*integral;
    else
        Xk(1+k) = (1/L)*integral;
    end
end

hold on;
Ntimes = 200;
t_eval = linspace(-L/2, L/2, Ntimes);
f_eval = f(t_eval);
plot(t_eval,f_eval)

%reconstrucci0n
fourier = zeros(1, Ntimes);
for k=ks
    onda = @(t)(exp(1j*(2*pi*k/L)*t));
    onda_eval = onda(t_eval);
    if k<0
        fourier = fourier + Xk(2*orden+2+k)*onda_eval;
    else
        fourier = fourier + Xk(1+k)*onda_eval;
    end
end

max(abs(imag(fourier)))%comprobamos que la parte imaginaria es cero
plot(t_eval, real(fourier), 'g-')
title('Reconstruccion con la serie de Fourier')
legends = cell(1, 2);
legends{1} = 'f';
legends{2} = sprintf('Serie de Fourier hasta grado %d',orden);
legend(legends)
hold off;

ans =    1.2439e-16

En esta otra gráfica, mostramos cómo la aproximación mejora al añadir más términos de la serie de Fourier.

In [5]:
f = @(x)(1 + exp(-x.^2).*cos(5.*x));

orden = 8;
ks = [-orden:orden];
Xk = zeros(2*orden+1);
for k = ks
    integrando = @(t) (f(t).*exp(-1j*(2*pi*k/L)*t));
    integral = complex_quadrature(integrando,-L/2,L/2);
    if k<0
        Xk(2*orden+2+k) = (1/L)*integral;
    else
        Xk(1+k) = (1/L)*integral;
    end
end

hold on;
legends = {};
Ntimes = 200;
t_eval = linspace(-L/2, L/2, Ntimes);
f_eval = f(t_eval);
plot(t_eval,f_eval)
lengends = cell(1, 1+orden);
legends{1} = 'f';

%reconstruccion
fourier = ones(1, Ntimes)*Xk(1);
for k=[1:orden]
    onda_plus = @(t)(exp(1j*(2*pi*k/L)*t));
    onda_plus_eval = onda_plus(t_eval);
    onda_minus = @(t)(exp(-1j*(2*pi*k/L)*t));
    onda_minus_eval = onda_minus(t_eval);
    fourier = fourier + Xk(2*orden+2-k)*onda_minus_eval + Xk(1+k)*onda_plus_eval;
    plot(t_eval, fourier)
    legends{k+1} = sprintf('Serie de Fourier hasta grado %d',k);
end

title('Reconstruccion con la serie de Fourier')
legend(legends);
hold off;

Fast Fourier Transform

El algoritmo conocido como Fast Fourier Transform permite calcular de forma eficiente simultáneamente todos los números: $$ X_k = \sum_{j=0}^{n-1} x_je^{-2\pi i k\frac{j}{n}} $$ dada una serie de números $x_j$.

La serie de números $\{X_k\}\;k=1\dots n$ se llama la transformada de Fourier discreta (en inglés Discrete Fourier Transform, o DFT) de la serie $\{x_k\}\;k=1\dots n$.

Vamos a usar la DFT para aproximar los coeficientes de la serie de Fourier de la función $x$, de esta forma:

  • Primero cambiamos de variable $t\in[0,L]$ a $\tau=\frac{t}{L}\in[0,1]$: $$ \hat{x}_k = \frac{1}{L}\int_{0}^L x(t)e^{-i2\pi \frac{k}{L} t}\:dt = \int_{0}^1 x(L\tau)e^{-i2\pi k \tau}\:d\tau $$
  • Después aproximamos usando el valor del integrando en $n$ puntos intermedios (es la regla del punto izquierdo, pero si la señal es periódica también admite otras interpretaciones, y es una aproximación mucho mejor) $$ \hat{f}_k \approx \frac{1}{n}\sum_{j=0}^{n-1} x(L \frac{j}{n}) e^{-i2\pi k \frac{j}{n}} $$
  • Observamos que si creamos un vector con $x_j=x(L \frac{j}{n})$, podemos obtener aproximaciones a todos los coeficientes de la serie de Fourier en una sóla llamada al algoritmo FFT: $$ \hat{x}_k \approx \frac{1}{n}X_k $$
In [6]:
L = 2*pi;
f = @(t)(1 + sign(sin(2*t)));

orden = 24;
ks = [-orden: orden];
h = L/(2*orden+1);
ts = linspace(0, L-h, 2*orden+1);
fs = f(ts);
Xk = fft(fs)*(1/(2*orden+1));

%print(Xk) %Imprime los coeficientes de la serie de Fourier
hold on;
ks_ordenado = [ks(orden+1:2*orden+1), ks(1:orden)];
plot(ks_ordenado, real(Xk),'go')
plot(ks_ordenado, imag(Xk),'bo')
legends = cell(1, 2);
legends{1} = 'parte real';
legends{2} =  'parte imaginaria';
legend(legends)
title('Coeficientes de la serie de Fourier')
hold off;
In [7]:
hold on;
Ntimes = 200;
t_eval = linspace(-L/2, L/2, Ntimes);
f_eval = f(t_eval);
plot(t_eval,f_eval)

%reconstruccion
fourier = zeros(1, Ntimes);
for k=ks
    onda = @(t)(exp(1j*(2*pi*k/L)*t));
    onda_eval = onda(t_eval);
    if k<0
        fourier = fourier + Xk(2*orden+2+k)*onda_eval;
    else
        fourier = fourier + Xk(1+k)*onda_eval;
    end
end

max(abs(imag(fourier)))%comprobamos que la parte imaginaria es cero
plot(t_eval, fourier, 'g-')
title('Reconstruccion con la serie de Fourier')
legend(sprintf('Serie de Fourier hasta grado %d',orden))
hold off;

ans =    3.2022e-16

Observamos que los coeficientes obtenidos mediante FFT son sólo aproximaciones a los auténticos coeficientes de la serie de Fourier.

In [8]:
L = 2*pi;
f = @(t)(1 + sign(sin(2*t)));

orden = 24;
% Usando quad
ks = [-orden: orden];
ks = [-orden:orden];
Xk = zeros(2*orden+1);
for k = ks
    integrando = @(t) (f(t).*exp(-1j*(2*pi*k/L)*t));
    integral = complex_quadrature(integrando,-L/2,L/2);
    if k<0
        Xk(2*orden+2+k) = (1/L)*integral;
    else
        Xk(1+k) = (1/L)*integral;
    end
end

% Usando fft
h = L/(2*orden+1);
ts = linspace(0, L-h, 2*orden+1);
fs = f(ts);
Xkf = fft(fs)*(1/(2*orden+1));

%print(Xk) %Imprime los coeficientes de la serie de Fourier
hold on;
%plot(ks, real(Xk),'go')#, label='parte real')
ks_ordenado = [ks(orden+1:2*orden+1), ks(1:orden)];
plot(ks_ordenado, real(Xk),'b.')
plot(ks_ordenado, imag(Xk),'k.')
plot(ks_ordenado, real(Xkf),'y+')
plot(ks_ordenado, imag(Xkf),'r+')
legends = cell(1, 4);
legends{1} = 'parte real';
legends{2} = 'parte imaginaria';
legends{3} = 'parte real (FFT)';
legends{4} = 'parte imaginaria (FFT)';
legend(legends)
title('Coeficientes de la serie de Fourier')
hold off;
In [9]:
hold on;
Ntimes = 200;
t_eval = linspace(-L/2, L/2, Ntimes);
f_eval = f(t_eval);
plot(t_eval,f_eval, 'r-')

%reconstruccion
fourier = zeros(1, Ntimes);
fourier_fft = zeros(1, Ntimes);
for k=ks
    onda = @(t)(exp(1j*(2*pi*k/L)*t));
    onda_eval = onda(t_eval);
    if k<0
        fourier = fourier + Xk(2*orden+2+k)*onda_eval;
        fourier_fft = fourier_fft + Xkf(2*orden+2+k)*onda_eval;
    else
        fourier = fourier + Xk(1+k)*onda_eval;
        fourier_fft = fourier_fft + Xkf(1+k)*onda_eval;
    end
end

plot(t_eval, fourier, 'g-')
plot(t_eval, fourier_fft, 'b-')
title('Reconstruccion con la serie de Fourier')
legends = cell(1, 3);
legends{1} = 'f';
legends{2} = sprintf('Serie de Fourier hasta grado %d',orden);
legends{3} = sprintf('Serie de Fourier hasta grado %d (FFT)',orden);
legend(legends)
hold off;

Sin embargo, los números obtenidos mediante FFT no se suelen interpretar como aproximaciones a los coeficientes de Fourier de una señal continua, sino como una transformación de una señal discreta. La transformada de Fourier tiene muchas aplicaciones en comunidades muy distintas, y es por eso que hay muchas variantes de las definiciones, a menudo con sutiles diferencias.

Transformada de Fourier de una señal no acotada

Recordamos que para una señal que no es periódica ni acotada, no se puede hablar de su serie de Fourier, pero sí de su transformada de Fourier (si la señal decae lo bastante rápido): $$ \hat{x}(w) = \int_{-\infty}^\infty x(t)e^{-i2\pi w t}\:dt $$ donde la reconstrucción es posible usando la fórmula de inversión: $$ x(t) = \int_{-\infty}^\infty \hat{x}(w)e^{i2\pi w t}\:dw $$

Ejemplo

Calculamos la transformada de Fourier y después reconstruimos una función a valores reales.

In [40]:
f = @(w)((1+10*w.^2).*exp( -(w-.1).^4));
L = 4;
xs = linspace(-2,2,200);
ss = f(xs);
plot(xs,ss)

Para calcular los valores de la transformada de Fourier, simplemente evaluamos las integrales de antes para unos cuantos valores equiespaciados de w: $w_k=k \: \Delta w$, para k desde -n hasta n.

Reemplazamos la integral desde $-\infty$ hasta $\infty$ por una integral finita (para eso es importante que la señal $x(t)$ decaiga rápidamente):

$$ \hat{x}(w) = \int_{-\infty}^\infty x(t)e^{-i2\pi w t}\:dt \approx \int_{-L/2}^{L/2} x(t)e^{-i2\pi w t}\:dt $$
In [57]:
n = 30;
Nw = 2*n;
wmax = 3;
ws = linspace(-wmax, wmax, Nw+1);

% Usando quad
Xk = zeros(1,Nw+1);
for k = [1:Nw+1]
    w = ws(k);
    integrando = @(t) (f(t).*exp(-1j*(2*pi*w)*t));
    integral = complex_quadrature(integrando,-L,L);
    Xk(k) = integral;
end

%print(Xk) %Imprime los coeficientes de la serie de Fourier
hold on;
plot(ws, real(Xk),'b')
plot(ws, imag(Xk),'k')
legends = {'parte real',
'parte imaginaria'};
hold off;
legend(legends)
title('Transformada de Fourier')

Ahora reconstruímos la función original de forma aproximada, usando una aproximación a la integral sobre w que aparece en la fórmula de inversión:

$$ x(t) = \int_{-\infty}^\infty \hat{x}(w)e^{i2\pi w t}\:dw \approx \int_{-w_{max}}^{w_{max}} \hat{x}(w)e^{i2\pi w t}\:dw \approx \Delta w \sum_{k=-n}^{n} \hat{x}(w_k)e^{i2\pi w_k t} $$

donde $w_{max} = n \Delta w $, y recorremos un total de $N=2n+1$ frecuencias: $$-n\Delta w,\dots,-\Delta w,0,\Delta w ,\dots,n\Delta w$$

In [56]:
hold on;
Ntimes = 200;
t_eval = linspace(-L/2, L/2, Ntimes);
f_eval = f(t_eval);
plot(t_eval,f_eval, 'r-')

%reconstruccion
dw = 2*wmax/Nw;
fourier = zeros(1, Ntimes);
for k=[1:Nw+1]
    w = ws(k);
    onda = @(t)(exp(1j*(2*pi*w)*t));
    onda_eval = onda(t_eval);
    fourier = fourier + dw*Xk(k)*onda_eval;
end

plot(t_eval, fourier, 'g.')
title('Reconstruccion con la Transformada de Fourier')
legends = cell(1, 2);
legends{1} = 'f';
legends{2} = 'Transformada de Fourier';
legend(legends)
hold off;

Los coeficientes se pueden aproximar usando la DFT:

  • Primero reemplazamos la integral desde $-\infty$ hasta $\infty$ por una integral finita, como hicimos antes (para eso es importante que decaiga rápidamente):
$$ \hat{x}(w) = \int_{-\infty}^\infty x(t)e^{-i2\pi w t}\:dt \approx \int_{-L/2}^{L/2} x(t)e^{-i2\pi w t}\:dt $$

  • Luego cambiamos de variable $t\in[-L/2,L/2]$ a $\tau=\frac{t}{L}+\frac{1}{2}\in[0,1]$, y aproximamos $w$ por $k/L$: $$ \hat{x}(w) \approx \hat{x}_k = \int_{-L/2}^{L/2} x(t)e^{-i2\pi \frac{k}{L} t}\:dt = L\int_{0}^1 x\left(L(\tau-1/2)\right)e^{-i2\pi k (\tau-1/2)}\:d\tau $$ $$ \hat{x}(w) \approx = L e^{i\pi k}\int_{0}^1 x\left(L(\tau-1/2)\right)e^{-i2\pi k \tau}\:d\tau $$

  • Después aproximamos usando el valor del integrando en $N$ puntos (esto hace un total de N-1 intervalos de igual longitud) $$ \hat{x}_k \approx L\frac{e^{i\pi k}}{N-1}\sum_{j=0}^{N-1} x\left(L \frac{j}{N} - \frac{L}{2}\right) e^{-i2\pi k \frac{j}{N}} $$

  • Observamos que si creamos un vector con $x_j=x(L \frac{j}{N}- \frac{L}{2})$, podemos obtener aproximaciones a los valores de $\hat{x}$ en todos los puntos de la forma $k/L$ en una sóla llamada al algoritmo FFT: $$ \hat{x}(k/L) \approx\hat{x}_k \approx \frac{Le^{i\pi k}}{N-1}X_k $$

Repetimos ahora el cálculo de la transformada de Fourier y la reconstrucción usando fft en vez de quad.

In [75]:
L = 30;
orden = 100;
% de este modo NW es impar
Nw = 2*orden + 1;

ts = linspace(-L/2, L/2, Nw);
fs = f(ts);
ks = [[0:orden],[-orden:-1]];
dw = 1/L;
Xkf = fft(fs).*exp(1j*pi*ks)*(L/(Nw-1));

%print(Xk) %Imprime los coeficientes de la serie de Fourier
hold on;
ws_ordenado = [-orden:orden]/L;
Xkf_ordenado = [Xkf(orden+1:Nw),Xkf(1:orden)];
plot(ws_ordenado, real(Xkf_ordenado),'g-')
plot(ws_ordenado, imag(Xkf_ordenado),'r-')
legends = {'parte real (FFT)',
'parte imaginaria (FFT)'};
hold off;
legend(legends)
title('Transformada de Fourier (FFT)')
In [83]:
hold on;
Ntimes = 200;
t_eval = linspace(-L/2, L/2, Ntimes);
f_eval = f(t_eval);
plot(t_eval,f_eval, 'r-')

%reconstruccion
dw = 1/L;
fourier_fft = zeros(1, Ntimes);
for k=[1:Nw]
    w = ws_ordenado(k);
    onda = @(t)(exp(1j*(2*pi*w)*t));
    onda_eval = onda(t_eval);
    fourier_fft = fourier_fft + dw*Xkf_ordenado(k)*onda_eval;
end

plot(t_eval, fourier_fft, 'k--')
title('Reconstruccion con la Transformada de Fourier')
legends = cell(1, 2);
legends{1} = 'f';
legends{2} = 'Transformada de Fourier';
legend(legends)
hold off;

Ejercicio

Usamos una función importante en el estudio del oleaje: la función de espectro de Pierson-Moskowitz. Aunque normalmente se define sólo para w>0, definimos también la versión definida para todo $w\in\mathbb{R}$, definida para w<0 mediante S(w)=S(-w), ya que la función de densidad espectral es una función par.

  • Calcula y dibuja la transformada de Fourier de espectro de Pierson-Moskowitz, en su versión definida en todos los reales.
  • Reconstruye la función para comprobar tu resultado.
  • Repite el cálculo usando fft y también integración con complex_quadrature.
In [84]:
%En esta versión del espectro de Pierson-Moskowitz 
% hay dos parámetros:
Hs = 5;  %Altura significativa de ola
Tp = 10; %Periodo pico de ola

%A partir de estos parametros se calculan A y B:
A = (5/16)*Hs^2*Tp^(-4);
B = (5/4)*Tp^(-4);

S = @(w)(A*w.^(-5) .* exp( -B*w.^(-4) ));

%Extendemos ahora el espectro de Pierson-Moskowitz
% a los numeros reales negativos por paridad
% Usamos la funcion de Heaviside
H = @(w) (w > 0);
S_bi = @(w)(H(w).*S(w) + H(-w).*S(-w));
%Ahora S_bi(-w) = S_bi(w)
In [85]:
ws = linspace(0.01,1,200);
ss = S(ws);
quad(S, 0,1)
plot(ws,ss);
title('Espectro de Pierson Moscowitz')

ans =  1.5623
In [86]:
ws = linspace(-1,1,200);
ss = S_bi(ws);
quad(S_bi, 0,1)
plot(ws,ss);
title('Espectro de Pierson Moscowitz')

ans =  1.5623
In [ ]:
In [ ]:
In [ ]:

Reconstrucción de Nyquist-Shannon

Recordamos que si tenemos una señal limitada en banda, la podemos reconstruir a partir de muestras equiespaciadas, si esas muestras se toman con suficiente frecuencia.

Si la señal es $x$, y tenemos muestras $x(hj)$, para $h$ igual al inverso de la frecuencia de muestreo, podemos reconstruir la señal de forma perfecta usando la fórmula de Nyquist-Shannon: $$ x(t)=\sum _{n=-\infty }^{\infty }x(nT)\cdot \mathrm {sinc} \left({\frac {t-nT}{T}}\right), $$ y la podemos recontruir de forma aproximada limitando la suma $$ x(t)\approx\sum _{n=-N }^{N }x(nT)\cdot \mathrm {sinc} \left({\frac {t-nT}{T}}\right), $$

In [15]:
%Dibujamos la funcion `sinc`
t_eval = linspace(-5, 5, 200);
f_eval = sinc(t_eval);
plot(t_eval,f_eval, 'DisplayName', 'f');
title('La funcion sinc')

Ejemplo

Esta señal es una suma de 4 ondas con frecuencias 1, 2, 3 y 6, con ciertas amplitudes, y donde cada onda se desfasa una cantidad aleatoria. De este modo, podemos seleccionar varias funciones diferentes, pero todas están acotadas en banda, con b=6.

In [16]:
function v=sumaondas(t, Fs, As, Phis)
    n = length(t);
    v = zeros(1, n);
    k = length(Fs);
    for j =[1:k]
        a = As(j);
        f = Fs(j);
        phi = Phis(j);
        v = v + a*cos(2*pi*f*t+ phi);
    end
end
In [17]:
T = 2;
Fs = [1, 2, 3, 6];
As = [1, 1.5, 1, 0.5];
Phis = 2*pi*rand(1,4);

% Limitada en ancho de banda
f = @(t)(sumaondas(t, Fs, As, Phis));

t_eval = linspace(0, T , 200);
f_eval = f(t_eval);
plot(t_eval,f_eval);
title('Onda aleatoria sumando ondas con fases aleatorias');

Intentamos reconstruir la función anterior, que era aleatoria, usando sólo unas muestras espaciadas $\frac{1}{2\cdot 6}$, que es la frecuencia de Nyquist. Observamos que un sampleado con mayor frecuencia que la frecuencia de Nyquist, la reconstrucción es perfecta, mientras que con un sampleado con menor frecuencia, la reconstrucción tiene fallos.

In [18]:
sample_interval = 1/12+0.01;
%sample_interval = 1/12-0.01
T = 2;
t_sample = [0:sample_interval:T];
Ntimes = 200;
t_eval = linspace(0,T, Ntimes);
f_eval = f(t_eval);
reconstruccion = zeros(1,Ntimes);

for k = [-floor(T/sample_interval):floor(2*T/sample_interval)]
    reconstruccion = reconstruccion + f(sample_interval*k)*sinc((t_eval - sample_interval*k)/sample_interval);
end

t_eval = linspace(0, T , 200);
f_eval = f(t_eval);
hold on
plot(t_eval,f_eval, 'b-')
plot(t_eval,reconstruccion, 'g')
hold off
legend('f', 'reconstruccion')

Oleaje generado con el espectro de Pierson-Moskowitz

Generamos ahora señales que son suma de varias ondas con frecuencias $0,dw,2dw,\cdots, 1$ equiespaciadas entre 0 y 1, con ciertas amplitudes que se derivan del espectro de Pierson-Moskowitz, y donde cada onda se desfasa una cantidad aleatoria.

De este modo, podemos seleccionar varias funciones diferentes, que corresponden a oleajes posibles, pero todas están acotadas en banda. En este caso el ancho de banda no está tan claramente definido, porque en principio hay ondas con cualquier frecuencia, aunque la amplitud que corresponde a frecuencias grandes es casi cero. De este modo, la frecuencia de Nyquist no está tan clara, pero puedes probar a cambiar la frecuencia y encontrarás rápidamente qué frecuencias permiten una reconstrucción (casi) perfecta. Observa también que al desfasar las ondas, no sólo tenemos frecuencias positivas sino también negativas.

In [19]:
T = 20;

n_ondas_simples = 20;
wepsilon = 1e-4;
wmax = 1;
dw = (wmax - wepsilon)/n_ondas_simples;
Fs = linspace(wepsilon,wmax, n_ondas_simples);
As = zeros(1, n_ondas_simples);
Phis = zeros(1, n_ondas_simples);
for j = [1:n_ondas_simples]
    w = Fs(j);
    As(j) = sqrt(2*S(w)*dw);
    Phis(j) = 2*pi*rand();
end

f = @(t)(sumaondas(t, Fs, As, Phis));

t_eval = linspace(0,T, 200);
f_eval = f(t_eval);

plot(t_eval , f_eval)
xlabel('time [ s ]')
ylabel('signal')
legend('a random wave')
In [20]:
T = 20;
sample_interval = 0.9;

t_sample = [0:sample_interval:T];
Ntimes = 200;
t_eval = linspace(0,T, Ntimes);
f_eval = f(t_eval);

reconstruccion = zeros(1, Ntimes);

for k = [-floor(T/sample_interval):floor(2*T/sample_interval)]
    reconstruccion = reconstruccion + f(sample_interval*k)*sinc((t_eval - sample_interval*k)/sample_interval);
end

hold on
plot(t_eval,f_eval, 'b-')
plot(t_eval,reconstruccion, 'g')
hold off
legend('input data', 'reconstruccion')
In [ ]: