Métodos no paramétricos
En las secciones siguientes se describen los métodos , , , y de estimación no paramétrica, junto con los , y .periodogramaperiododegrafía modificadaWelchmultitaperFunción CPSDestimación de la 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 transformación de Fourier en tiempo discreto de las muestras del proceso (normalmente se realiza en una cuadrícula con un FFT) y escalar adecuadamente la magnitud cuadrada del resultado. Esta estimación se denomina .periodograma
La estimación del periodo de la psD de una señal
donde Fs es la frecuencia de muestreo.
En la práctica, el cómputo real de
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 1001 elementos, 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 de periodograma de la PSD se puede calcular utilizando .periodograma En este caso, el vector de datos se multiplica por una ventana 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')
Algoritmo
El periodograma calcula y escala la salida de la FFT para producir la gráfica de potencia frente a la frecuencia de la siguiente manera.
Si la señal de entrada tiene un valor real, la magnitud del FFT resultante es simétrica con respecto a la frecuencia cero (DC). Para un FFT de longitud uniforme, solo los primeros (1 + /2) puntos son únicos.
nfftDetermine el número de valores únicos y mantenga solo esos puntos únicos.Tome las magnitudes cuadradas de los valores FFT únicos. Escalar las magnitudes cuadradas (excepto para DC) por
<math display="block">, donde está la longitud de la señal antes de cualquier relleno cero.N Escalar el valor de CC en<mrow></math><mn>2</mn><mo>/</mo><mo stretchy="false">(</mo><msub><mrow><mi>F</mi></mrow><mrow></msub><mi>s</mi></mrow><mi>N</mi><mo stretchy="false">)</mo></mrow><math display="block">.<mrow></math><mn>1</mn><mo>/</mo><mo stretchy="false">(</mo><msub><mrow><mi>F</mi></mrow><mrow></msub><mi>s</mi></mrow><mi>N</mi><mo stretchy="false">)</mo></mrow>Cree un vector de frecuencia a partir del número de puntos únicos, el nfft y la frecuencia de muestreo.
Trazar la magnitud resultante al cuadrado FFT contra la frecuencia.
Rendimiento del Periodograma
En las secciones siguientes se analiza el rendimiento del periodograma con respecto a los problemas de fugas, resolución, sesgo y varianza.
Fuga espectral
Considere el PSD de una longitud finita (longitud
Dado 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
mostrando que el valor esperado del periodograma es la convolución del verdadero PSD con el cuadrado del núcleo Dirichlet.
El efecto de la convolución se entiende mejor para los datos sinusoidales. Supongamos que
Su espectro es
que para una secuencia de longitud finita se convierte en
La ecuación anterior es igual a
Así que en el espectro de la señal de longitud finita, los deltas de Dirac han sido reemplazados por términos de la forma
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 parcela muestra un lóbulo principal y varios lóbulos 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 potencia concentrada exactamente en los puntos de frecuencia discretos
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 de 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 depende únicamente de 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.
Resolución
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.Resolución
Para resolver dos sinusoides que están relativamente cerca entre sí en 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 estos sinusoides. La anchura del lóbulo principal se define como la anchura 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). Este ancho es aproximadamente igual a
En otras palabras, para dos sinusoides de frecuencias
En el ejemplo anterior, donde dos sinusoides están separados por solo 10 Hz, el registro de datos debe ser superior a 100 muestras para permitir la resolución de dos sinusoides distintos por agram.
Considere un caso en el que no se cumpla este criterio, en cuanto a la secuencia de 67 muestras a continuación:
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 ahora. Cuando el SNR es bajo, las características espectrales verdaderas son mucho más difíciles de distinguir, y los artefactos de ruido aparecen en estimaciones espectrales basadas en el periodograma. El ejemplo siguiente 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)
Sesgo del Periodograma
El periodograma es un estimador sesgado del PSD. Anteriormente se demostró que su valor esperado era
El periodograma es asintóticamente imparcial, lo que es evidente desde 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 a la función delta de Dirac. Sin embargo, en algunos casos el periodo de periodo según el periodo es un estimador deficiente de la PSD, incluso cuando el registro de datos es largo. Esto se debe a la varianza del periodograma, como se explica a continuación.
Variación del Periodograma
Se puede demostrar que la varianza del periodograma
lo que indica que la varianza no tiende a cero como la longitud de los datos
El Periodograma Modificado
Las ventanas la señal del dominio del tiempo antes de calcular el DFT para suavizar los bordes de la señal.periododegrafía modificada Esto tiene el efecto de reducir la altura de los lóbulos laterales o fugas espectrales. Este fenómeno da lugar a la interpretación de los lóbulos laterales como frecuencias espurias 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 graves. Por otro lado, las ventanas no rectangulares también amplían el lóbulo principal, lo que se traduce en una reducción de la resolución.
Le 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 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)
Se puede comprobar que aunque los lóbulos laterales son mucho menos evidentes en el periodograma de ventana Hamming, los dos picos principales son más anchos. De hecho, el ancho de 3 dB del lóbulo principal correspondiente a una ventana hamming es aproximadamente el doble que el de una ventana rectangular. Por lo tanto, para una longitud de datos fija, la resolución PSD alcanzable con una ventana Hamming es aproximadamente la mitad que se puede lograr con una ventana rectangular. Los intereses de la competencia de la anchura del lóbulo principal y la altura del lóbulo lateral se pueden resolver en cierta medida mediante el uso de ventanas variables como la ventana Kaiser.
Las ventanas no rectangulares afectan a la potencia media de una señal porque algunas de las muestras de tiempo se atenúan cuando se multiplican por la ventana. Para compensar esto, y normalizar la ventana para tener un poder promedio de unidad.periodogramapwelch Esto garantiza que la potencia media medida sea generalmente independiente de la elección de la ventana. Si los componentes de frecuencia no son bien resueltos por los estimadores PSD, la opción de ventana afecta a la potencia media.
La estimación de periodograma modificado de la PSD es
donde está la constante de normalización de la ventana:U
Para valores grandes de , tiende a ser independiente de la longitud de la ventana.LU La adición de como constante de normalización garantiza que el periodograma modificado sea asintóticamente imparcial.U
Método de Welch
Un estimador mejorado de la PSD es el propuesto por Welch. El método consiste en dividir los datos de la 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 de caja de herramientas implementa el método de Welch.pwelch
El promedio de los periodogramas modificados tiende a disminuir la varianza de la estimación en relación con una estimación de un solo periodograma de todo el registro de datos. Aunque la superposición entre segmentos introduce información redundante, este efecto se ve disminuido por el uso de una ventana no rectangular, lo que reduce la importancia o se da a las muestras finales de segmentos (las muestras que se superponen).Peso
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 varianza y la resolución. Se pueden 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 50% de superposición usando una ventana rectangular.
pwelch(xn,rectwin(150),50,512,fs)
En el periodo de arriba, el ruido y la fuga hacen que uno de los sinusoides sea esencialmente indistinguible de los picos artificiales. Por el contrario, aunque el PSD producido por el método de Welch tiene picos más amplios, todavía se pueden distinguir los dos sinusoides, que se destacan del "suelo de ruido".
Sin embargo, si tratamos de reducir aún más la varianza, la pérdida de resolución hace que uno de los 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 de la PSD. El valor esperado de la estimación de PSD es:
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 transformación de Fourier de la función de ventana.LUW(f) Como es el caso de todos los periodogramas, el estimador de Welch es asintó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 segmentos. Básicamente, la varianza es inversamente proporcional al número de segmentos cuyos periodogramas modificados se están promediando.
Método Multitaper
El periododeo se puede interpretar como filtro de una longitud
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 de PSD constante sobre el ancho de banda del filtro. Esto proporciona otra interpretación de por qué la estimación PSD del periodograma mejora a medida que aumenta la longitud de la señal. Sin embargo, hay dos factores evidentes desde 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álculo 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 puede recibir 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. De este modo, el modelo de banco de filtros proporciona una nueva interpretación del compromiso entre la varianza y la resolución.
Thompson's (MTM) se basa en estos resultados para proporcionar una estimación mejorada de PSD.método multitaper En lugar de utilizar filtros de paso de banda que son esencialmente ventanas rectangulares (como en el método de 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 esferoidales de prolato discretos (DPSS, también conocidos como ).Secuencias de Slepian
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 lo proporciona el producto de ancho de banda,
La función de ™ cuadro de herramientas de procesamiento de señales 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 esferoidales de prolato discretos. Para series de datos largas (10.000 puntos o más), es útil calcular los DPSS una vez y guardarlos en un archivo MAT. , , , y se proporcionan para mantener una base de datos de DSS guardados en el archivo MAT .dpsssavedpssloaddpssdirdpsscleardpss.mat
Función de densidad multiespectral
El PSD es un caso especial de la función (CPSD), definida entre dos señales ( ) y ( ) comodensidad espectral cruzadaxnyn
Como es el caso de las secuencias de correlación y covarianza, la caja de herramientas PSD y CPSD porque las longitudes de señal son finitas.Estimaciones
Para estimar la densidad transversal de dos señales de igual longitud y utilizando el método de Welch, la función forma el periodograma como el producto de la FFT y el conjugado de la FFT de .xycpsdxy A diferencia del PSD de valor real, el CPSD es una función compleja. maneja el seccionamiento y las ventanas de y de la misma manera que la función:cpsdxypwelch
Sxy = cpsd(x,y,nwin,noverlap,nfft,fs)
Estimación de la función de transferencia
Una aplicación del método de Welch es la identificación no paramétrica del sistema. Supongamos que se trata de un sistema lineal, invariable en el 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
Una estimación de la función de transferencia entre ( ) y ( ) esxnyn
Este método estima la información de magnitud y fase. La función utiliza el método de Welch para calcular el CPSD y el espectro de potencia, y luego forma su cociente para la estimación de la función de transferencia.tfestimate Utilice de la misma manera que utiliza la función.tfestimatecpsd
Generar una señal que consiste en dos sinusoides incrustados en el ruido gaussiano blanco.
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));
Filtrar la señal con un filtro de media móvil FIR.xn Calcular 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);Trazar 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
Este cociente es un número real entre 0 y 1 que mide la correlación entre ( ) y ( ) a la frecuenciaxnyn
La función toma secuencias y , calcula sus espectros de potencia y CPSD, y devuelve el cociente de la magnitud cuadrada de la CPSD y el producto de los espectros de potencia.mscoherexnyn Sus opciones y funcionamiento son similares a las funciones.cpsdtfestimate
Generar una señal que consiste en dos sinusoides incrustados en el ruido gaussiano blanco. 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));
Filtrar la señal con un filtro de media móvil FIR.xn Calcular y trazar 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 que.mscohere Esto se debe a que la función de coherencia para los datos dependientes linealmente es una.
Consulte también
Apps
Funciones
cpsd|mscohere|periodograma|pmtm|pwelch|tfestimate
Select a Web Site
Choose a web site to get translated content where available and see local events and offers. Based on your location, we recommend that you select: .
Select web siteYou can also select a web site from the following list:
Americas
- América Latina (Español)
- Canada (English)
- United States (English)
Europe
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)
