Esta página aún no se ha traducido para esta versión. Puede ver la versión más reciente de esta página en inglés.

Métodos no paramétricos

En las secciones siguientes se analizan los métodos de estimación no paramétrico, y, junto con los relacionados, y.periodogramaperiodograma modificadoWelchmulticonicidadFunción CPSDestimación de función de transferenciafunción de coherencia

Periodograma

En términos generales, una forma de estimar el PSD de un proceso es simplemente encontrar la transformada de Fourier en tiempo discreto de las muestras del proceso (usualmente hecho en una rejilla con una FFT) y escalar apropiadamente la magnitud cuadrada del resultado. Esta estimación se llama la.periodograma

La estimación del periodograma del PSD de una señal

<math display="block">
<mrow>
<msub>
<mrow>
<mi>x</mi>
</mrow>
<mrow>
<mi>L</mi>
</mrow>
</msub>
<mo stretchy="false">(</mo>
<mi>n</mi>
<mo stretchy="false">)</mo>
</mrow>
</math>
de longitud esL

<math display="block">
<mrow>
<msub>
<mrow>
<mi>P</mi>
</mrow>
<mrow>
<mi>x</mi>
<mi>x</mi>
</mrow>
</msub>
<mo stretchy="false">(</mo>
<mi>f</mi>
<mo stretchy="false">)</mo>
<mo>=</mo>
<mrow>
<mfrac>
<mrow>
<mn>1</mn>
</mrow>
<mrow>
<mrow>
<mi>L</mi>
<msub>
<mrow>
<mi>F</mi>
</mrow>
<mrow>
<mi>s</mi>
</mrow>
</msub>
</mrow>
</mrow>
</mfrac>
</mrow>
<msup>
<mrow>
<mo>|</mo>
<munderover>
<mrow>
<mo></mo>
</mrow>
<mrow>
<mi>n</mi>
<mo>=</mo>
<mn>0</mn>
</mrow>
<mrow>
<mi>L</mi>
<mo>-</mo>
<mn>1</mn>
</mrow>
</munderover>
<msub>
<mrow>
<mi>x</mi>
</mrow>
<mrow>
<mi>L</mi>
</mrow>
</msub>
<mo stretchy="false">(</mo>
<mi>n</mi>
<mo stretchy="false">)</mo>
<msup>
<mrow>
<mi>e</mi>
</mrow>
<mrow>
<mo>-</mo>
<mi>j</mi>
<mn>2</mn>
<mi>π</mi>
<mi>f</mi>
<mi>n</mi>
<mo>/</mo>
<msub>
<mrow>
<mi>F</mi>
</mrow>
<mrow>
<mi>s</mi>
</mrow>
</msub>
</mrow>
</msup>
<mo>|</mo>
</mrow>
<mrow>
<mn>2</mn>
</mrow>
</msup>
<mo>,</mo>
</mrow>
</math>

donde FS es la frecuencia de muestreo.

En la práctica, el cómputo real de

<math display="block">
<mrow>
<msub>
<mrow>
<mi>P</mi>
</mrow>
<mrow>
<mi>x</mi>
<mi>x</mi>
</mrow>
</msub>
<mo stretchy="false">(</mo>
<mi>f</mi>
<mo stretchy="false">)</mo>
</mrow>
</math>
sólo se puede realizar en un número finito de puntos de frecuencia, y por lo general emplea una FFT. La mayoría de las implementaciones del método de periodograma computan el
<math display="block">
<mrow>
<mi>N</mi>
</mrow>
</math>
estimación de PSD de punto en las frecuencias

<math display="block">
<mrow>
<msub>
<mrow>
<mi>f</mi>
</mrow>
<mrow>
<mi>k</mi>
</mrow>
</msub>
<mo>=</mo>
<mrow>
<mfrac>
<mrow>
<mrow>
<mi>k</mi>
<msub>
<mrow>
<mi>F</mi>
</mrow>
<mrow>
<mi>s</mi>
</mrow>
</msub>
</mrow>
</mrow>
<mrow>
<mi>N</mi>
</mrow>
</mfrac>
</mrow>
<mo>,</mo>
<mspace width="0.2777777777777778em"></mspace>
<mspace width="0.2777777777777778em"></mspace>
<mspace width="0.2777777777777778em"></mspace>
<mi>k</mi>
<mo>=</mo>
<mn>0</mn>
<mo>,</mo>
<mn>1</mn>
<mo>,</mo>
<mo></mo>
<mo>,</mo>
<mi>N</mi>
<mo>-</mo>
<mn>1</mn>
<mo>.</mo>
</mrow>
</math>

En algunos casos, el cálculo del periodograma a través de un algoritmo FFT es más eficiente si el número de frecuencias es una potencia de dos. Por lo tanto, no es raro rellenar la señal de entrada con ceros para extender su longitud a una potencia de dos.

Como ejemplo del periodograma, considere la siguiente señal de elemento 1001, que consta de dos sinusoides más ruido:xn

fs = 1000;                % Sampling frequency t = (0:fs)/fs;            % One second worth of samples A = [1 2];                % Sinusoid amplitudes (row vector) f = [150;140];            % Sinusoid frequencies (column vector) xn = A*sin(2*pi*f*t) + 0.1*randn(size(t)); % The three last lines are equivalent to % xn = sin(2*pi*150*t) + 2*sin(2*pi*140*t) + 0.1*randn(size(t));

La estimación del periodograma del PSD se puede calcular utilizando.periodograma En este caso, el vector de datos se multiplica por una ventana de Hamming para producir un periodograma modificado.

[Pxx,F] = periodogram(xn,hamming(length(xn)),length(xn),fs); plot(F,10*log10(Pxx)) xlabel('Hz') ylabel('dB') title('Modified Periodogram Power Spectral Density Estimate')

Algorithm

Periodograma calcula y escala la salida de la FFT para producir el diagrama de potencia frente a la frecuencia de la siguiente manera.

  1. Si la señal de entrada es de valor real, la magnitud de la FFT resultante es simétrica con respecto a la frecuencia cero (DC). Para una FFT de longitud incluso, solo los primeros puntos (1 +/2) son únicos.nfft Determine el número de valores únicos y conserve solo esos puntos únicos.

  2. Tome las magnitudes cuadradas de los valores FFT únicos. Escale las magnitudes cuadradas (excepto para DC) por

    <math display="block">
    <mrow>
    <mn>2</mn>
    <mo>/</mo>
    <mo stretchy="false">(</mo>
    <msub>
    <mrow>
    <mi>F</mi>
    </mrow>
    <mrow>
    <mi>s</mi>
    </mrow>
    </msub>
    <mi>N</mi>
    <mo stretchy="false">)</mo>
    </mrow>
    </math>
    , donde es la longitud de la señal antes de cualquier relleno de cero.N Escale el valor de DC
    <math display="block">
    <mrow>
    <mn>1</mn>
    <mo>/</mo>
    <mo stretchy="false">(</mo>
    <msub>
    <mrow>
    <mi>F</mi>
    </mrow>
    <mrow>
    <mi>s</mi>
    </mrow>
    </msub>
    <mi>N</mi>
    <mo stretchy="false">)</mo>
    </mrow>
    </math>
    .

  3. Cree un vector de frecuencia a partir del número de puntos únicos, el nfft y la frecuencia de muestreo.

  4. Graficar la magnitud resultante al cuadrado FFT contra la frecuencia.

El rendimiento del periodograma

En las secciones siguientes se analiza el rendimiento del periodograma con respecto a los problemas de fuga, resolución, sesgo y varianza.

Spectral Leakage

Considere el PSD de una longitud finita (longitud

<math display="block">
<mrow>
<mi>L</mi>
</mrow>
</math>
) señal
<math display="block">
<mrow>
<msub>
<mrow>
<mi>x</mi>
</mrow>
<mrow>
<mi>L</mi>
</mrow>
</msub>
<mo stretchy="false">(</mo>
<mi>n</mi>
<mo stretchy="false">)</mo>
</mrow>
</math>
. Con frecuencia es útil interpretar
<math display="block">
<mrow>
<msub>
<mrow>
<mi>x</mi>
</mrow>
<mrow>
<mi>L</mi>
</mrow>
</msub>
<mo stretchy="false">(</mo>
<mi>n</mi>
<mo stretchy="false">)</mo>
</mrow>
</math>
como resultado de multiplicar una señal infinita,
<math display="block">
<mrow>
<mi>x</mi>
<mo stretchy="false">(</mo>
<mi>n</mi>
<mo stretchy="false">)</mo>
</mrow>
</math>
, por una ventana rectangular de longitud finita,
<math display="block">
<mrow>
<msub>
<mrow>
<mi>w</mi>
</mrow>
<mrow>
<mi>R</mi>
</mrow>
</msub>
<mo stretchy="false">(</mo>
<mi>n</mi>
<mo stretchy="false">)</mo>
</mrow>
</math>
:

<math display="block">
<mrow>
<msub>
<mrow>
<mi>x</mi>
</mrow>
<mrow>
<mi>L</mi>
</mrow>
</msub>
<mo stretchy="false">(</mo>
<mi>n</mi>
<mo stretchy="false">)</mo>
<mo>=</mo>
<mi>x</mi>
<mo stretchy="false">(</mo>
<mi>n</mi>
<mo stretchy="false">)</mo>
<msub>
<mrow>
<mi>w</mi>
</mrow>
<mrow>
<mi>R</mi>
</mrow>
</msub>
<mo stretchy="false">(</mo>
<mi>n</mi>
<mo stretchy="false">)</mo>
<mo>.</mo>
</mrow>
</math>

Debido a que la multiplicación en el dominio de tiempo corresponde a la convolución en el dominio de frecuencia, el valor esperado del periodograma en el dominio de frecuencia es

<math display="block">
<mrow>
<mi>E</mi>
<mo stretchy="false">{</mo>
<mrow>
<msub>
<mrow>
<munderover accent="true">
<mrow>
<mi>P</mi>
</mrow>
<mrow></mrow>
<mrow>
<mo>ˆ</mo>
</mrow>
</munderover>
</mrow>
<mrow>
<mi>x</mi>
<mi>x</mi>
</mrow>
</msub>
</mrow>
<mo stretchy="false">(</mo>
<mi>f</mi>
<mo stretchy="false">)</mo>
<mo stretchy="false">}</mo>
<mo>=</mo>
<mrow>
<mfrac>
<mrow>
<mn>1</mn>
</mrow>
<mrow>
<mrow>
<mrow>
<msub>
<mrow>
<mi>F</mi>
</mrow>
<mrow>
<mi>s</mi>
</mrow>
</msub>
</mrow>
</mrow>
</mrow>
</mfrac>
</mrow>
<msubsup>
<mrow>
<mo></mo>
</mrow>
<mrow>
<mo>-</mo>
<mrow>
<msub>
<mrow>
<mi>F</mi>
</mrow>
<mrow>
<mi>s</mi>
</mrow>
</msub>
</mrow>
<mo>/</mo>
<mn>2</mn>
</mrow>
<mrow>
<mrow>
<msub>
<mrow>
<mi>F</mi>
</mrow>
<mrow>
<mi>s</mi>
</mrow>
</msub>
</mrow>
<mo>/</mo>
<mn>2</mn>
</mrow>
</msubsup>
<mrow>
<mrow>
<mfrac>
<mrow>
<mrow>
<mrow>
<msup>
<mrow>
<mi mathvariant="normal">sin</mi>
</mrow>
<mrow>
<mn>2</mn>
</mrow>
</msup>
</mrow>
<mo stretchy="false">(</mo>
<mi>L</mi>
<mi>π</mi>
<mo stretchy="false">(</mo>
<mi>f</mi>
<mo>-</mo>
<msup>
<mrow>
<mi>f</mi>
</mrow>
<mrow>
<mo></mo>
</mrow>
</msup>
<mo stretchy="false">)</mo>
<mo>/</mo>
<mrow>
<msub>
<mrow>
<mi>F</mi>
</mrow>
<mrow>
<mi>s</mi>
</mrow>
</msub>
</mrow>
<mo stretchy="false">)</mo>
</mrow>
</mrow>
<mrow>
<mrow>
<mi>L</mi>
<mrow>
<msup>
<mrow>
<mi mathvariant="normal">sin</mi>
</mrow>
<mrow>
<mn>2</mn>
</mrow>
</msup>
</mrow>
<mo stretchy="false">(</mo>
<mi>π</mi>
<mo stretchy="false">(</mo>
<mi>f</mi>
<mo>-</mo>
<msup>
<mrow>
<mi>f</mi>
</mrow>
<mrow>
<mo></mo>
</mrow>
</msup>
<mo stretchy="false">)</mo>
<mo>/</mo>
<mrow>
<msub>
<mrow>
<mi>F</mi>
</mrow>
<mrow>
<mi>s</mi>
</mrow>
</msub>
</mrow>
<mo stretchy="false">)</mo>
</mrow>
</mrow>
</mfrac>
</mrow>
</mrow>
<mspace width="0.16666666666666666em"></mspace>
<mrow>
<msub>
<mrow>
<mi>P</mi>
</mrow>
<mrow>
<mi>x</mi>
<mi>x</mi>
</mrow>
</msub>
</mrow>
<mo stretchy="false">(</mo>
<msup>
<mrow>
<mi>f</mi>
</mrow>
<mrow>
<mo></mo>
</mrow>
</msup>
<mo stretchy="false">)</mo>
<mspace width="0.16666666666666666em"></mspace>
<mi>d</mi>
<msup>
<mrow>
<mi>f</mi>
</mrow>
<mrow>
<mo></mo>
</mrow>
</msup>
<mo>,</mo>
</mrow>
</math>

mostrando que el valor esperado del periodograma es la convolución del PSD verdadero con el cuadrado del kernel de Dirichlet.

El efecto de la convolución se entiende mejor para los datos sinusoidales. Supongamos que

<math display="block">
<mrow>
<mi>x</mi>
<mo stretchy="false">(</mo>
<mi>n</mi>
<mo stretchy="false">)</mo>
</mrow>
</math>
se compone de una suma de
<math display="block">
<mrow>
<mi>M</mi>
</mrow>
</math>
sinusoides complejos:

<math display="block">
<mrow>
<mi>x</mi>
<mo stretchy="false">(</mo>
<mi>n</mi>
<mo stretchy="false">)</mo>
<mo>=</mo>
<munderover>
<mrow>
<mo></mo>
</mrow>
<mrow>
<mi>k</mi>
<mo>=</mo>
<mn>1</mn>
</mrow>
<mrow>
<mi>N</mi>
</mrow>
</munderover>
<mrow>
<mrow>
<msub>
<mrow>
<mi>A</mi>
</mrow>
<mrow>
<mi>k</mi>
</mrow>
</msub>
</mrow>
</mrow>
<mrow>
<msup>
<mrow>
<mi>e</mi>
</mrow>
<mrow>
<mi>j</mi>
<mrow>
<msub>
<mrow>
<mi>ω</mi>
</mrow>
<mrow>
<mi>k</mi>
</mrow>
</msub>
</mrow>
<mi>n</mi>
</mrow>
</msup>
</mrow>
<mo>.</mo>
</mrow>
</math>

Su espectro es

<math display="block">
<mrow>
<mi>X</mi>
<mo stretchy="false">(</mo>
<mi>ω</mi>
<mo stretchy="false">)</mo>
<mo>=</mo>
<munderover>
<mrow>
<mo></mo>
</mrow>
<mrow>
<mi>k</mi>
<mo>=</mo>
<mn>1</mn>
</mrow>
<mrow>
<mi>N</mi>
</mrow>
</munderover>
<mrow>
<msub>
<mrow>
<mi>A</mi>
</mrow>
<mrow>
<mi>k</mi>
</mrow>
</msub>
</mrow>
<mi>δ</mi>
<mo stretchy="false">(</mo>
<mi>ω</mi>
<mo>-</mo>
<mrow>
<msub>
<mrow>
<mi>ω</mi>
</mrow>
<mrow>
<mi>k</mi>
</mrow>
</msub>
</mrow>
<mo stretchy="false">)</mo>
<mo>,</mo>
</mrow>
</math>

que para una secuencia de longitud finita se convierte en

<math display="block">
<mrow>
<mi>X</mi>
<mo stretchy="false">(</mo>
<mi>ω</mi>
<mo stretchy="false">)</mo>
<mo>=</mo>
<msubsup>
<mrow>
<mo></mo>
</mrow>
<mrow>
<mo>-</mo>
<mi>π</mi>
</mrow>
<mrow>
<mi>π</mi>
</mrow>
</msubsup>
<mrow>
<mrow>
<mo>(</mo>
<mrow>
<munderover>
<mrow>
<mo></mo>
</mrow>
<mrow>
<mi>k</mi>
<mo>=</mo>
<mn>1</mn>
</mrow>
<mrow>
<mi>N</mi>
</mrow>
</munderover>
<mrow>
<mrow>
<msub>
<mrow>
<mi>A</mi>
</mrow>
<mrow>
<mi>k</mi>
</mrow>
</msub>
</mrow>
<mi>δ</mi>
<mo stretchy="false">(</mo>
<mi>ε</mi>
<mo>-</mo>
<mrow>
<msub>
<mrow>
<mi>ω</mi>
</mrow>
<mrow>
<mi>k</mi>
</mrow>
</msub>
</mrow>
<mo stretchy="false">)</mo>
</mrow>
</mrow>
<mo>)</mo>
</mrow>
<mrow>
<msub>
<mrow>
<mi>W</mi>
</mrow>
<mrow>
<mi>R</mi>
</mrow>
</msub>
</mrow>
<mo stretchy="false">(</mo>
<mi>ω</mi>
<mo>-</mo>
<mi>ε</mi>
<mo stretchy="false">)</mo>
<mspace width="0.16666666666666666em"></mspace>
<mi>d</mi>
<mi>ε</mi>
</mrow>
<mo>.</mo>
</mrow>
</math>

La ecuación anterior es igual a

<math display="block">
<mrow>
<mi>X</mi>
<mo stretchy="false">(</mo>
<mi>ω</mi>
<mo stretchy="false">)</mo>
<mo>=</mo>
<munderover>
<mrow>
<mo></mo>
</mrow>
<mrow>
<mi>k</mi>
<mo>=</mo>
<mn>1</mn>
</mrow>
<mrow>
<mi>N</mi>
</mrow>
</munderover>
<mrow>
<mrow>
<msub>
<mrow>
<mi>A</mi>
</mrow>
<mrow>
<mi>k</mi>
</mrow>
</msub>
</mrow>
</mrow>
<mrow>
<msub>
<mrow>
<mi>W</mi>
</mrow>
<mrow>
<mi>R</mi>
</mrow>
</msub>
</mrow>
<mo stretchy="false">(</mo>
<mi>ω</mi>
<mo>-</mo>
<mrow>
<msub>
<mrow>
<mi>ω</mi>
</mrow>
<mrow>
<mi>k</mi>
</mrow>
</msub>
</mrow>
<mo stretchy="false">)</mo>
<mo>.</mo>
</mrow>
</math>

Por lo tanto, en el espectro de la señal de longitud finita, los deltas Dirac han sido reemplazados por términos de la forma

<math display="block">
<mrow>
<msub>
<mrow>
<mi>W</mi>
</mrow>
<mrow>
<mi>R</mi>
</mrow>
</msub>
<mo stretchy="false">(</mo>
<mi>ω</mi>
<mo>-</mo>
<msub>
<mrow>
<mi>ω</mi>
</mrow>
<mrow>
<mi>k</mi>
</mrow>
</msub>
<mo stretchy="false">)</mo>
</mrow>
</math>
, que corresponde a la respuesta de frecuencia de una ventana rectangular centrada en la frecuencia
<math display="block">
<mrow>
<msub>
<mrow>
<mi>ω</mi>
</mrow>
<mrow>
<mi>k</mi>
</mrow>
</msub>
</mrow>
</math>
.

La respuesta de frecuencia de una ventana rectangular tiene la forma de un sinc periódico:

L = 32; [h,w] = freqz(rectwin(L)/L,1); y = diric(w,L);  plot(w/pi,20*log10(abs(h))) hold on plot(w/pi,20*log10(abs(y)),'--') hold off ylim([-40,0]) legend('Frequency Response','Periodic Sinc') xlabel('\omega / \pi')

La trama muestra un lóbulo principal y varios laterales, el más grande de los cuales es aproximadamente 13,5 dB por debajo del pico del lóbulo principal. Estos lóbulos representan el efecto conocido como fuga espectral. Mientras que la señal de longitud infinita tiene su poder concentrado exactamente en los puntos de frecuencia discretos

<math display="block">
<mrow>
<msub>
<mrow>
<mi>f</mi>
</mrow>
<mrow>
<mi>k</mi>
</mrow>
</msub>
</mrow>
</math>
, la señal de ventana (o truncada) tiene un continuo de energía "filtrada" alrededor de los puntos de frecuencia discretos
<math display="block">
<mrow>
<msub>
<mrow>
<mi>f</mi>
</mrow>
<mrow>
<mi>k</mi>
</mrow>
</msub>
</mrow>
</math>
.

Debido a que la respuesta de frecuencia de una ventana rectangular corta es una aproximación mucho más pobre a la función Delta Dirac que la de una ventana más larga, la fuga espectral es especialmente evidente cuando los registros de datos son cortos. Considere la siguiente secuencia de 100 muestras:

fs = 1000;                 % Sampling frequency t = (0:fs/10)/fs;          % One-tenth second worth of samples A = [1 2];                 % Sinusoid amplitudes f = [150;140];             % Sinusoid frequencies xn = A*sin(2*pi*f*t) + 0.1*randn(size(t)); periodogram(xn,rectwin(length(xn)),1024,fs)

Es importante tener en cuenta que el efecto de la fuga espectral está supeditado únicamente a la longitud del registro de datos. No es una consecuencia del hecho de que el periodograma se calcula en un número finito de muestras de frecuencia.

Resolution

se refiere a la capacidad de discriminar las características espectrales, y es un concepto clave en el análisis del rendimiento del estimador espectral.Resolution

Para resolver dos sinusoides que están relativamente cerca juntos en la frecuencia, es necesario que la diferencia entre las dos frecuencias sea mayor que la anchura del lóbulo principal de los espectros filtrados para cualquiera de estas sinusoides. El ancho del lóbulo principal se define como el ancho del lóbulo principal en el punto donde la potencia es la mitad de la potencia máxima del lóbulo principal (es decir, 3 dB de ancho). Esta anchura es aproximadamente igual a

<math display="block">
<mrow>
<msub>
<mrow>
<mi>f</mi>
</mrow>
<mrow>
<mi>s</mi>
</mrow>
</msub>
<mo>/</mo>
<mi>L</mi>
</mrow>
</math>
.

En otras palabras, para dos sinusoides de frecuencias

<math display="block">
<mrow>
<msub>
<mrow>
<mi>f</mi>
</mrow>
<mrow>
<mn>1</mn>
</mrow>
</msub>
</mrow>
</math>
Y
<math display="block">
<mrow>
<msub>
<mrow>
<mi>f</mi>
</mrow>
<mrow>
<mn>2</mn>
</mrow>
</msub>
</mrow>
</math>
, la condición de resolvabilidad requiere que

<math display="block">
<mrow>
<msub>
<mrow>
<mi>f</mi>
</mrow>
<mrow>
<mn>2</mn>
</mrow>
</msub>
<mo>-</mo>
<msub>
<mrow>
<mi>f</mi>
</mrow>
<mrow>
<mn>1</mn>
</mrow>
</msub>
<mo>></mo>
<mrow>
<mfrac>
<mrow>
<mrow>
<msub>
<mrow>
<mi>F</mi>
</mrow>
<mrow>
<mi>s</mi>
</mrow>
</msub>
</mrow>
</mrow>
<mrow>
<mi>L</mi>
</mrow>
</mfrac>
</mrow>
<mo>.</mo>
</mrow>
</math>

En el ejemplo anterior, donde dos sinusoides están separados por sólo 10 Hz, el registro de datos debe ser mayor que 100 muestras para permitir la resolución de dos sinusoides distintos por un periodograma.

Considere un caso en el que este criterio no se cumple, como para la secuencia de 67 muestras siguientes:

fs = 1000;                  % Sampling frequency t = (0:fs/15)/fs;           % 67 samples A = [1 2];                  % Sinusoid amplitudes f = [150;140];              % Sinusoid frequencies xn = A*sin(2*pi*f*t) + 0.1*randn(size(t)); periodogram(xn,rectwin(length(xn)),1024,fs)

La discusión anterior sobre la resolución no tuvo en cuenta los efectos del ruido, ya que la relación señal/ruido (SNR) ha sido relativamente alta hasta el momento. Cuando el SNR es bajo, las características espectrales verdaderas son mucho más difíciles de distinguir, y los artefactos del ruido aparecen en las estimaciones espectrales basadas en el periodograma. El siguiente ejemplo ilustra esto:

fs = 1000;                  % Sampling frequency t = (0:fs/10)/fs;           % One-tenth second worth of samples A = [1 2];                  % Sinusoid amplitudes f = [150;140];              % Sinusoid frequencies xn = A*sin(2*pi*f*t) + 2*randn(size(t)); periodogram(xn,rectwin(length(xn)),1024,fs)

Bias of the Periodogram

El periodograma es un estimador sesgado del PSD. Su valor esperado fue mostrado anteriormente como

<math display="block">
<mrow>
<mi>E</mi>
<mo stretchy="false">{</mo>
<mrow>
<msub>
<mrow>
<munderover accent="true">
<mrow>
<mi>P</mi>
</mrow>
<mrow></mrow>
<mrow>
<mo>ˆ</mo>
</mrow>
</munderover>
</mrow>
<mrow>
<mi>x</mi>
<mi>x</mi>
</mrow>
</msub>
</mrow>
<mo stretchy="false">(</mo>
<mi>f</mi>
<mo stretchy="false">)</mo>
<mo stretchy="false">}</mo>
<mo>=</mo>
<mrow>
<mfrac>
<mrow>
<mn>1</mn>
</mrow>
<mrow>
<mrow>
<mrow>
<msub>
<mrow>
<mi>F</mi>
</mrow>
<mrow>
<mi>s</mi>
</mrow>
</msub>
</mrow>
</mrow>
</mrow>
</mfrac>
</mrow>
<msubsup>
<mrow>
<mo></mo>
</mrow>
<mrow>
<mo>-</mo>
<mrow>
<msub>
<mrow>
<mi>F</mi>
</mrow>
<mrow>
<mi>s</mi>
</mrow>
</msub>
</mrow>
<mo>/</mo>
<mn>2</mn>
</mrow>
<mrow>
<mrow>
<msub>
<mrow>
<mi>F</mi>
</mrow>
<mrow>
<mi>s</mi>
</mrow>
</msub>
</mrow>
<mo>/</mo>
<mn>2</mn>
</mrow>
</msubsup>
<mrow>
<mrow>
<mfrac>
<mrow>
<mrow>
<mrow>
<msup>
<mrow>
<mi mathvariant="normal">sin</mi>
</mrow>
<mrow>
<mn>2</mn>
</mrow>
</msup>
</mrow>
<mo stretchy="false">(</mo>
<mi>L</mi>
<mi>π</mi>
<mo stretchy="false">(</mo>
<mi>f</mi>
<mo>-</mo>
<msup>
<mrow>
<mi>f</mi>
</mrow>
<mrow>
<mo></mo>
</mrow>
</msup>
<mo stretchy="false">)</mo>
<mo>/</mo>
<mrow>
<msub>
<mrow>
<mi>F</mi>
</mrow>
<mrow>
<mi>s</mi>
</mrow>
</msub>
</mrow>
<mo stretchy="false">)</mo>
</mrow>
</mrow>
<mrow>
<mrow>
<mi>L</mi>
<mrow>
<msup>
<mrow>
<mi mathvariant="normal">sin</mi>
</mrow>
<mrow>
<mn>2</mn>
</mrow>
</msup>
</mrow>
<mo stretchy="false">(</mo>
<mi>π</mi>
<mo stretchy="false">(</mo>
<mi>f</mi>
<mo>-</mo>
<msup>
<mrow>
<mi>f</mi>
</mrow>
<mrow>
<mo></mo>
</mrow>
</msup>
<mo stretchy="false">)</mo>
<mo>/</mo>
<mrow>
<msub>
<mrow>
<mi>F</mi>
</mrow>
<mrow>
<mi>s</mi>
</mrow>
</msub>
</mrow>
<mo stretchy="false">)</mo>
</mrow>
</mrow>
</mfrac>
</mrow>
</mrow>
<mspace width="0.16666666666666666em"></mspace>
<mrow>
<msub>
<mrow>
<mi>P</mi>
</mrow>
<mrow>
<mi>x</mi>
<mi>x</mi>
</mrow>
</msub>
</mrow>
<mo stretchy="false">(</mo>
<msup>
<mrow>
<mi>f</mi>
</mrow>
<mrow>
<mo></mo>
</mrow>
</msup>
<mo stretchy="false">)</mo>
<mspace width="0.16666666666666666em"></mspace>
<mi>d</mi>
<msup>
<mrow>
<mi>f</mi>
</mrow>
<mrow>
<mo></mo>
</mrow>
</msup>
<mo>.</mo>
</mrow>
</math>

El periodograma es asintóticamente imparcial, lo cual es evidente a partir de la observación anterior que a medida que la longitud del registro de datos tiende al infinito, la respuesta de frecuencia de la ventana rectangular se aproxima más estrechamente a la función Delta del Dirac. Sin embargo, en algunos casos el periodograma es un pobre estimador del PSD incluso cuando el registro de datos es largo. Esto se debe a la varianza del periodograma, como se explica a continuación.

Variance of the Periodogram

Se puede demostrar que la varianza del periodograma es

<math display="block">
<mrow>
<mrow>
<mstyle mathvariant="normal">
<mrow>
<mi>V</mi>
<mi>a</mi>
<mi>r</mi>
</mrow>
</mstyle>
</mrow>
<mo stretchy="false">(</mo>
<msub>
<mrow>
<munderover accent="true">
<mrow>
<mi>P</mi>
</mrow>
<mrow></mrow>
<mrow>
<mo>ˆ</mo>
</mrow>
</munderover>
</mrow>
<mrow>
<mi>x</mi>
<mi>x</mi>
</mrow>
</msub>
<mo stretchy="false">(</mo>
<mi>f</mi>
<mo stretchy="false">)</mo>
<mo stretchy="false">)</mo>
<mo>=</mo>
<mrow>
<mo>{</mo>
<mtable>
<mtr>
<mtd>
<mrow>
<msubsup>
<mrow>
<mi>P</mi>
</mrow>
<mrow>
<mi>x</mi>
<mi>x</mi>
</mrow>
<mrow>
<mn>2</mn>
</mrow>
</msubsup>
<mo stretchy="false">(</mo>
<mi>f</mi>
<mo stretchy="false">)</mo>
<mo>,</mo>
</mrow>
</mtd>
<mtd>
<mrow>
<mn>0</mn>
<mo><</mo>
<mi>f</mi>
<mo><</mo>
<msub>
<mrow>
<mi>F</mi>
</mrow>
<mrow>
<mi>s</mi>
</mrow>
</msub>
<mo>/</mo>
<mn>2</mn>
<mo>,</mo>
</mrow>
</mtd>
</mtr>
<mtr>
<mtd>
<mrow>
<mn>2</mn>
<msubsup>
<mrow>
<mi>P</mi>
</mrow>
<mrow>
<mi>x</mi>
<mi>x</mi>
</mrow>
<mrow>
<mn>2</mn>
</mrow>
</msubsup>
<mo stretchy="false">(</mo>
<mi>f</mi>
<mo stretchy="false">)</mo>
<mo>,</mo>
</mrow>
</mtd>
<mtd>
<mrow>
<mi>f</mi>
<mo>=</mo>
<mn>0</mn>
<mo>,</mo>
<mspace width="0.2777777777777778em"></mspace>
<msub>
<mrow>
<mi>F</mi>
</mrow>
<mrow>
<mi>s</mi>
</mrow>
</msub>
<mo>/</mo>
<mn>2</mn>
<mo>,</mo>
</mrow>
</mtd>
</mtr>
</mtable>
</mrow>
</mrow>
</math>

que indica que la varianza no tiende a cero como la longitud de datos

<math display="block">
<mrow>
<mi>L</mi>
</mrow>
</math>
tiende al infinito. En términos estadísticos, el periodograma no es un estimador consistente del PSD. Sin embargo, el periodograma puede ser una herramienta útil para la estimación espectral en situaciones donde el SNR es alto, y especialmente si el registro de datos es largo.

El periodograma modificado

Las ventanas de la señal de dominio de tiempo antes de calcular el DFT con el fin de suavizar los bordes de la señal.periodograma modificado Esto tiene el efecto de reducir la altura de la margen lateral o fuga espectral. Este fenómeno da lugar a la interpretación de las bandas laterales como frecuencias no esenciales introducidas en la señal por el truncamiento abrupto que se produce cuando se utiliza una ventana rectangular. Para las ventanas no rectangulares, los puntos finales de la señal truncada se atenúan suavemente, y por lo tanto las frecuencias no esenciales introducidas son mucho menos severas. Por otro lado, las ventanas no rectangulares también amplían el lóbulo principal, lo que da como resultado una reducción de la resolución.

El permite calcular un periodograma modificado especificando la ventana que se utilizará en los datos.periodograma Por ejemplo, compare una ventana rectangular predeterminada y una ventana de Hamming. Especifique el mismo número de puntos DFT en ambos casos.

fs = 1000;                  % Sampling frequency t = (0:fs/10)/fs;           % One-tenth second worth of samples A = [1 2];                  % Sinusoid amplitudes f = [150;140];              % Sinusoid frequencies nfft = 1024;  xn = A*sin(2*pi*f*t) + 0.1*randn(size(t)); periodogram(xn,rectwin(length(xn)),nfft,fs)

periodogram(xn,hamming(length(xn)),nfft,fs)

Usted puede verificar que aunque los laterales son mucho menos evidentes en el periodograma Hamming-windowed, los dos picos principales son más anchos. De hecho, la anchura de 3 dB del lóbulo principal correspondiente a una ventana de Hamming es aproximadamente el doble que la de una ventana rectangular. Por lo tanto, para una longitud de datos fija, la resolución PSD alcanzable con una ventana de Hamming es aproximadamente la mitad que alcanzable con una ventana rectangular. Los intereses competidores de ancho del lóbulo principal y altura de mínimos se pueden resolver hasta cierto punto mediante el uso de ventanas de variables como la ventana de Kaiser.

Las ventanas no rectangulares afectan a la potencia media de una señal porque algunas de las muestras de tiempo se atenen cuando se multiplican por la ventana. Para compensar esto, y normalizar la ventana para tener una potencia media de la unidad.periodogramapwelch Esto garantiza que la potencia media medida es generalmente independiente de la opción de ventana. Si los componentes de frecuencia no están bien resueltos por los estimadores PSD, la opción de ventana afecta a la potencia media.

La estimación del periodograma modificado del PSD es

<math display="block">
<mrow>
<msub>
<mrow>
<munderover accent="true">
<mrow>
<mi>P</mi>
</mrow>
<mrow></mrow>
<mrow>
<mo>ˆ</mo>
</mrow>
</munderover>
</mrow>
<mrow>
<mi>x</mi>
<mi>x</mi>
</mrow>
</msub>
<mo stretchy="false">(</mo>
<mi>f</mi>
<mo stretchy="false">)</mo>
<mo>=</mo>
<mrow>
<mfrac>
<mrow>
<mrow>
<msup>
<mrow>
<mo>|</mo>
<mi>X</mi>
<mo stretchy="false">(</mo>
<mi>f</mi>
<mo stretchy="false">)</mo>
<mo>|</mo>
</mrow>
<mrow>
<mn>2</mn>
</mrow>
</msup>
</mrow>
</mrow>
<mrow>
<mrow>
<msub>
<mrow>
<mi>F</mi>
</mrow>
<mrow>
<mi>s</mi>
</mrow>
</msub>
<mi>L</mi>
<mi>U</mi>
</mrow>
</mrow>
</mfrac>
</mrow>
<mo>,</mo>
</mrow>
</math>

donde está la constante de normalización de ventana:U

<math display="block">
<mrow>
<mi>U</mi>
<mo>=</mo>
<mrow>
<mfrac>
<mrow>
<mn>1</mn>
</mrow>
<mrow>
<mi>L</mi>
</mrow>
</mfrac>
</mrow>
<munderover>
<mrow>
<mo></mo>
</mrow>
<mrow>
<mi>n</mi>
<mo>=</mo>
<mn>0</mn>
</mrow>
<mrow>
<mi>N</mi>
<mo>-</mo>
<mn>1</mn>
</mrow>
</munderover>
<msup>
<mrow>
<mo>|</mo>
<mi>w</mi>
<mo stretchy="false">(</mo>
<mi>n</mi>
<mo stretchy="false">)</mo>
<mo>|</mo>
</mrow>
<mrow>
<mn>2</mn>
</mrow>
</msup>
<mo>.</mo>
</mrow>
</math>

Para los valores grandes de, tiende a ser independiente de la longitud de la ventana.LU La adición de una constante de normalización asegura que el periodograma modificado sea asintóticamente imparcial.U

El método de Welch

Un estimador mejorado del PSD es el propuesto por Welch. El método consiste en dividir los datos de serie temporal en segmentos (posiblemente superpuestos), calcular un periodograma modificado de cada segmento y, a continuación, promediar las estimaciones de PSD. El resultado es la estimación de PSD de Welch. La función Toolbox implementa el método de Welch.pwelch

El promedio de periodogramas modificados tiende a disminuir la varianza de la estimación en relación con una sola estimación de periodograma de todo el registro de datos. Aunque la superposición entre segmentos introduce información redundante, este efecto se disminuye por el uso de una ventana no rectangular, lo que reduce la importancia o se da a las muestras finales de los segmentos (los ejemplos que se superponen).weight

Sin embargo, como se mencionó anteriormente, el uso combinado de registros de datos cortos y ventanas no rectangulares da como resultado una resolución reducida del estimador. En Resumen, hay un equilibrio entre la reducción de la varianza y la resolución. Uno puede manipular los parámetros en el método de Welch para obtener estimaciones mejoradas en relación con el periodograma, especialmente cuando el SNR es bajo. Esto se ilustra en el ejemplo siguiente.

Considere una señal que consta de 301 muestras:

fs = 1000;             % Sampling frequency t = (0:0.3*fs)/fs;     % 301 samples A = [2 8];             % Sinusoid amplitudes (row vector) f = [150;140];         % Sinusoid frequencies (column vector)  xn = A*sin(2*pi*f*t) + 5*randn(size(t)); periodogram(xn,rectwin(length(xn)),1024,fs)

Podemos obtener la estimación espectral de Welch para 3 segmentos con un 50% de superposición usando una ventana rectangular.

pwelch(xn,rectwin(150),50,512,fs)

En el periodograma anterior, el ruido y la fuga hacen que una de las sinusoides sea esencialmente indistinguible de los picos artificiales. En contraste, aunque el PSD producido por el método de Welch tiene picos más anchos, todavía se pueden distinguir las dos sinusoides, que destacan desde el "piso de ruido."

Sin embargo, si tratamos de reducir aún más la varianza, la pérdida de resolución hace que una de las sinusoides se pierda por completo.

pwelch(xn,rectwin(100),75,512,fs)

Sesgo y normalización en el método de Welch

El método de Welch produce un estimador sesgado del PSD. El valor esperado de la estimación PSD es:

E{PWelch(f)}=1FsLUFs/2Fs/2|W(ff)|2Pxx(f)df,

donde está la longitud de los segmentos de datos, es la misma constante de normalización presente en la definición del periodograma modificado, y es la transformada de Fourier de la función de ventana.LUW(f) Como es el caso de todos los periodogramas, el estimador de Welch es asintomáticamente imparcial. Para un registro de datos de longitud fija, el sesgo de la estimación de Welch es mayor que el del periodograma porque la longitud de los segmentos es menor que la longitud de toda la muestra de datos.

La varianza del estimador de Welch es difícil de calcular porque depende tanto de la ventana utilizada como de la cantidad de superposición entre los segmentos. Básicamente, la varianza es inversamente proporcional al número de segmentos cuyos periodogramas modificados se están promediando.

Método multiconicidad

El periodograma puede interpretarse como un filtrado de una longitud

<math display="block">
<mrow>
<mi>L</mi>
</mrow>
</math>
Señal
<math display="block">
<mrow>
<msub>
<mrow>
<mi>x</mi>
</mrow>
<mrow>
<mi>L</mi>
</mrow>
</msub>
<mo stretchy="false">(</mo>
<mi>n</mi>
<mo stretchy="false">)</mo>
</mrow>
</math>
, a través de un banco filtrante (un conjunto de filtros en paralelo) de
<math display="block">
<mrow>
<mi>L</mi>
</mrow>
</math>
Filtros de paso de banda FIR. El ancho de banda de 3 dB de cada uno de estos filtros de paso de banda se puede demostrar que es aproximadamente igual a
<math display="block">
<mrow>
<msub>
<mrow>
<mi>f</mi>
</mrow>
<mrow>
<mi>s</mi>
</mrow>
</msub>
<mo>/</mo>
<mi>L</mi>
</mrow>
</math>
. La respuesta de magnitud de cada uno de estos filtros de paso de banda se asemeja a la de una ventana rectangular. Por lo tanto, el periodograma se puede ver como un cálculo de la potencia de cada señal filtrada (es decir, la salida de cada filtro de paso de banda) que utiliza sólo una muestra de cada señal filtrada y asume que el PSD de
<math display="block">
<mrow>
<msub>
<mrow>
<mi>x</mi>
</mrow>
<mrow>
<mi>L</mi>
</mrow>
</msub>
<mo stretchy="false">(</mo>
<mi>n</mi>
<mo stretchy="false">)</mo>
</mrow>
</math>
es constante sobre el ancho de banda de cada filtro de paso de banda.

A medida que aumenta la longitud de la señal, el ancho de banda de cada filtro de paso de banda disminuye, lo que lo convierte en un filtro más selectivo y mejora la aproximación del PSD constante sobre el ancho de banda del filtro. Esto proporciona otra interpretación de por qué la estimación de PSD del periodograma mejora a medida que aumenta la longitud de la señal. Sin embargo, hay dos factores que se desprenden de este punto de vista que comprometen la precisión de la estimación del periodograma. En primer lugar, la ventana rectangular produce un filtro de paso de banda deficiente. En segundo lugar, el cómputo de la potencia en la salida de cada filtro de paso de banda se basa en una sola muestra de la señal de salida, produciendo una aproximación muy cruda.

El método de Welch se puede dar una interpretación similar en términos de un banco de filtros. En la implementación de Welch, se utilizan varias muestras para calcular la potencia de salida, lo que resulta en una variación reducida de la estimación. Por otro lado, el ancho de banda de cada filtro de paso de banda es mayor que el correspondiente al método de periodograma, lo que resulta en una pérdida de resolución. El modelo de banco de filtros proporciona así una nueva interpretación del compromiso entre la varianza y la resolución.

Thompson (MTM) se basa en estos resultados para proporcionar una estimación de PSD mejorada.multitaper method En lugar de utilizar filtros de paso de banda que son esencialmente ventanas rectangulares (como en el método periodograma), el método MTM utiliza un banco de filtros de paso de banda óptimos para calcular la estimación. Estos filtros FIR óptimos se derivan de un conjunto de secuencias conocidas como secuencias esferidales prolato discretas (DPSSs, también conocido como).Slepian sequences

Además, el método MTM proporciona un parámetro de ancho de banda de tiempo con el que equilibrar la varianza y la resolución. Este parámetro viene dado por el producto de ancho de banda de tiempo,

<math display="block">
<mrow>
<mi>N</mi>
<mi>W</mi>
</mrow>
</math>
y está directamente relacionado con el número de se estrecha utilizados para computar el espectro. Siempre hay
<math display="block">
<mrow>
<mn>2</mn>
<mi>N</mi>
<mi>W</mi>
<mo>-</mo>
<mn>1</mn>
</mrow>
</math>
se estrecha utilizados para formar la estimación. Esto significa que, como
<math display="block">
<mrow>
<mi>N</mi>
<mi>W</mi>
</mrow>
</math>
aumenta, hay más estimaciones del espectro de potencia, y la varianza de la estimación disminuye. Sin embargo, el ancho de banda de cada conicidad también es proporcional a
<math display="block">
<mrow>
<mi>N</mi>
<mi>W</mi>
</mrow>
</math>
, así como
<math display="block">
<mrow>
<mi>N</mi>
<mi>W</mi>
</mrow>
</math>
aumenta, cada estimación exhibe más fugas espectrales (es decir, picos más anchos) y la estimación espectral general es más sesgada. Para cada conjunto de datos, normalmente hay un valor para
<math display="block">
<mrow>
<mi>N</mi>
<mi>W</mi>
</mrow>
</math>
que permite un óptimo intercambio entre sesgo y varianza.

La función™ de la caja de herramientas de procesamiento de señal que implementa el método MTM es.pmtm Se utiliza para calcular el PSD de una señal.pmtm

fs = 1000;                % Sampling frequency t = (0:fs)/fs;            % One second worth of samples A = [1 2];                % Sinusoid amplitudes f = [150;140];            % Sinusoid frequencies  xn = A*sin(2*pi*f*t) + 0.1*randn(size(t)); pmtm(xn,4,[],fs)

Al reducir el producto de ancho de banda de tiempo, puede aumentar la resolución a expensas de una varianza mayor.

pmtm(xn,1.5,[],fs)

Este método es más costoso computacionalmente que el método de Welch debido al costo de calcular las secuencias esferidales prolato discretas. Para series de datos largas (10.000 puntos o más), es útil calcular los DPSSs una vez y guardarlos en un archivo MAT. , y se proporcionan para mantener una base de datos de DPSSs guardados en el archivo MAT.dpsssavedpssloaddpssdirdpsscleardpss.mat

Función de densidad espectral cruzada

El PSD es un caso especial de la función (CPSD), definida entre dos señales () y () comocross spectral densityxnyn

Pxy(ω)=12πm=Rxy(m)ejωm.

Como es el caso de las secuencias de correlación y covarianza, el cuadro de herramientas el PSD y CPSD porque las longitudes de la señal son finitas.estimates

Para estimar la densidad espectral transversal de dos señales de igual longitud y utilizar el método de Welch, la función forma el periodograma como el producto de la FFT de y el conjugada de la FFT de.xycpsdxy A diferencia del PSD de valor real, el CPSD es una función compleja. maneja el corte y la ventana de la misma manera que la función:cpsdxypwelch

Sxy = cpsd(x,y,nwin,noverlap,nfft,fs) 

Función de transferencia estimación

Una aplicación del método de Welch es la identificación del sistema no paramétrico. Supongamos que es un sistema lineal, invariable de tiempo, y () y () son la entrada y salida de, respectivamente.HxnynH A continuación, el espectro de potencia de () está relacionado con el CPSD de () y () porxnxnyn

<math display="block">
<mrow>
<msub>
<mrow>
<mi>P</mi>
</mrow>
<mrow>
<mi>y</mi>
<mi>x</mi>
</mrow>
</msub>
<mo stretchy="false">(</mo>
<mi>ω</mi>
<mo stretchy="false">)</mo>
<mo>=</mo>
<mi>H</mi>
<mo stretchy="false">(</mo>
<mi>ω</mi>
<mo stretchy="false">)</mo>
<msub>
<mrow>
<mi>P</mi>
</mrow>
<mrow>
<mi>x</mi>
<mi>x</mi>
</mrow>
</msub>
<mo stretchy="false">(</mo>
<mi>ω</mi>
<mo stretchy="false">)</mo>
<mo>.</mo>
</mrow>
</math>

Una estimación de la función de transferencia entre () y () esxnyn

<math display="block">
<mrow>
<munderover accent="true">
<mrow>
<mi>H</mi>
</mrow>
<mrow></mrow>
<mrow>
<mo>ˆ</mo>
</mrow>
</munderover>
<mo stretchy="false">(</mo>
<mi>ω</mi>
<mo stretchy="false">)</mo>
<mo>=</mo>
<mrow>
<mfrac>
<mrow>
<mrow>
<msub>
<mrow>
<munderover accent="true">
<mrow>
<mi>P</mi>
</mrow>
<mrow></mrow>
<mrow>
<mo>ˆ</mo>
</mrow>
</munderover>
</mrow>
<mrow>
<mi>y</mi>
<mi>x</mi>
</mrow>
</msub>
<mo stretchy="false">(</mo>
<mi>ω</mi>
<mo stretchy="false">)</mo>
</mrow>
</mrow>
<mrow>
<mrow>
<msub>
<mrow>
<munderover accent="true">
<mrow>
<mi>P</mi>
</mrow>
<mrow></mrow>
<mrow>
<mo>ˆ</mo>
</mrow>
</munderover>
</mrow>
<mrow>
<mi>x</mi>
<mi>x</mi>
</mrow>
</msub>
<mo stretchy="false">(</mo>
<mi>ω</mi>
<mo stretchy="false">)</mo>
</mrow>
</mrow>
</mfrac>
</mrow>
<mo>.</mo>
</mrow>
</math>

Este método estima tanto la magnitud como la información de fase. La función utiliza el método de Welch para computar el CPSD y el espectro de potencia, y luego forma su cociente para la estimación de la función de transferencia.tfestimate Utilice la misma forma en que utiliza la función.tfestimatecpsd

Generar una señal consistente en dos sinusoides incrustados en el ruido blanco Gaussiano.

rng('default')  fs = 1000;                % Sampling frequency t = (0:fs)/fs;            % One second worth of samples A = [1 2];                % Sinusoid amplitudes f = [150;140];            % Sinusoid frequencies  xn = A*sin(2*pi*f*t) + 0.1*randn(size(t));

Filtre la señal con un filtro de media móvil FIR.xn Calcule la respuesta de magnitud real y la respuesta estimada.

h = ones(1,10)/10;             % Moving-average filter yn = filter(h,1,xn);  [HEST,f] = tfestimate(xn,yn,256,128,256,fs); H = freqz(h,1,f,fs);

Graficar los resultados.

subplot(2,1,1) plot(f,abs(H)) title('Actual Transfer Function Magnitude') yl = ylim; grid subplot(2,1,2) plot(f,abs(HEST)) title('Transfer Function Magnitude Estimate') xlabel('Frequency (Hz)') ylim(yl) grid

Función de coherencia

La coherencia de magnitud cuadrada entre dos señales () y () esxnyn

<math display="block">
<mrow>
<msub>
<mrow>
<mi>C</mi>
</mrow>
<mrow>
<mi>x</mi>
<mi>y</mi>
</mrow>
</msub>
<mo stretchy="false">(</mo>
<mi>ω</mi>
<mo stretchy="false">)</mo>
<mo>=</mo>
<mrow>
<mfrac>
<mrow>
<mrow>
<msup>
<mrow>
<mo>|</mo>
<msub>
<mrow>
<mi>P</mi>
</mrow>
<mrow>
<mi>x</mi>
<mi>y</mi>
</mrow>
</msub>
<mo stretchy="false">(</mo>
<mi>ω</mi>
<mo stretchy="false">)</mo>
<mo>|</mo>
</mrow>
<mrow>
<mn>2</mn>
</mrow>
</msup>
</mrow>
</mrow>
<mrow>
<mrow>
<msub>
<mrow>
<mi>P</mi>
</mrow>
<mrow>
<mi>x</mi>
<mi>x</mi>
</mrow>
</msub>
<mo stretchy="false">(</mo>
<mi>ω</mi>
<mo stretchy="false">)</mo>
<msub>
<mrow>
<mi>P</mi>
</mrow>
<mrow>
<mi>y</mi>
<mi>y</mi>
</mrow>
</msub>
<mo stretchy="false">(</mo>
<mi>ω</mi>
<mo stretchy="false">)</mo>
</mrow>
</mrow>
</mfrac>
</mrow>
<mo>.</mo>
</mrow>
</math>

Este cociente es un número real entre 0 y 1 que mide la correlación entre () y () en la frecuenciaxnyn

<math display="block">
<mrow>
<mi>ω</mi>
</mrow>
</math>
.

La función toma secuencias y, computa sus espectros de potencia y CPSD, y devuelve el cociente de la magnitud cuadrada del CPSD y el producto de los espectros de potencia.mscoherexnyn Sus opciones y funcionamiento son similares a las y funciones.cpsdtfestimate

Generar una señal consistente en dos sinusoides incrustados en el ruido blanco Gaussiano. La señal se muestrea a 1 kHz durante 1 segundo.

rng('default')  fs = 1000; t = (0:fs)/fs; A = [1 2];                % Sinusoid amplitudes f = [150;140];            % Sinusoid frequencies  xn = A*sin(2*pi*f*t) + 0.1*randn(size(t));

Filtre la señal con un filtro de media móvil FIR.xn Calcule y trace la función de coherencia y la salida del filtro en función de la frecuencia.xnyn

h = ones(1,10)/10; yn = filter(h,1,xn);  mscohere(xn,yn,256,128,256,fs)

Si la longitud de la secuencia de entrada, la longitud de la ventana y el número de puntos de datos superpuestos en una ventana son tales que funcionan en un solo registro, la función devuelve todos los.mscohere Esto se debe a que la función de coherencia para los datos dependientes linealmente es una.

Consulte también

Aplicaciones

Funciones