cpsd

Densidad espectral de potencia cruzada

Sintaxis

pxy = cpsd(x,y)

pxy = cpsd(x,y,window)

pxy = cpsd(x,y,window,noverlap)

pxy = cpsd(x,y,window,noverlap,nfft)

pxy = cpsd(___,'mimo')

[pxy,w] = cpsd(___)

[pxy,f] = cpsd(___,fs)

[pxy,w] = cpsd(x,y,window,noverlap,w)

[pxy,f] = cpsd(x,y,window,noverlap,f,fs)

[___] = cpsd(x,y,___,freqrange)

cpsd(___)

Descripción

pxy = cpsd(x,y) estima la densidad espectral de potencia cruzada (CPSD) de dos señales de tiempo discreto, x e y, utilizando el método de estimación espectral de periodograma modificado y promediado de Welch.

• Si x e y son ambos vectores, deben tener la misma longitud.

• Si una de las señales es una matriz y la otra es un vector, la longitud del vector debe ser igual al número de filas de la matriz. La función expande el vector y devuelve una matriz de estimaciones de densidad espectral de potencia cruzada columna por columna.

• Si x e y son matrices con el mismo número de filas, pero diferente número de columnas, cpsd devuelve un arreglo tridimensional, pxy, que contiene estimaciones de densidad espectral de potencia cruzada para todas las combinaciones de columnas de entrada. Cada columna de pxy corresponde a una columna de x, y cada página corresponde a una columna de y: pxy(:,m,n) = cpsd(x(:,m),y(:,n)).

• Si x e y son matrices del mismo tamaño, cpsd opera por columnas: pxy(:,n) = cpsd(x(:,n),y(:,n)). Para obtener un arreglo multi-entrada/multi-salida, añada 'mimo' a la lista de argumentos.

Para x e y reales, cpsd devuelve una CPSD unilateral. Para x o y complejos, cpsd devuelve una CPSD de dos caras.

pxy = cpsd(x,y,window) utiliza window para dividir x e y en segmentos y realizar ventaneo.

pxy = cpsd(x,y,window,noverlap) utiliza muestras de solapamiento noverlap entre segmentos contiguos.

pxy = cpsd(x,y,window,noverlap,nfft) utiliza puntos de muestreo nfft para calcular la transformada discreta de Fourier.

pxy = cpsd(___,'mimo') calcula un arreglo multi-entrada/multi-salida de estimaciones de densidad espectral de potencia cruzada. Esta sintaxis puede incluir cualquier combinación de argumentos de entrada de las sintaxis anteriores.

[pxy,w] = cpsd(___) devuelve un vector de frecuencias normalizadas, w, en las que se estima la densidad espectral de potencia cruzada.

[pxy,f] = cpsd(___,fs) devuelve un vector de frecuencias, f, expresado en términos de la tasa de muestreo, fs, a la que se estima la densidad espectral de potencia cruzada. fs debe ser la sexta entrada numérica a cpsd. Para introducir una tasa de muestreo y seguir utilizando los valores predeterminados de los argumentos opcionales anteriores, especifique estos argumentos como vacíos, [].

[pxy,w] = cpsd(x,y,window,noverlap,w) devuelve las estimaciones de densidad espectral de potencia cruzada en las frecuencias normalizadas especificadas en w.

[pxy,f] = cpsd(x,y,window,noverlap,f,fs) devuelve las estimaciones de densidad espectral de potencia cruzada en las frecuencias especificadas en f.

[___] = cpsd(x,y,___,freqrange) devuelve la estimación de la densidad espectral de potencia cruzada en el rango de frecuencias especificado por freqrange. Las opciones válidas para freqrange son 'onesided', 'twosided' y 'centered'.

cpsd(___) sin argumentos de salida representa la estimación de la densidad espectral de potencia cruzada en la ventana de la figura actual.

Ejemplos

Densidad espectral de potencia cruzada de señales de ruido de color

Genere dos señales de ruido de color y represente su densidad espectral de potencia cruzada. Especifique una FFT de longitud 1024 y una ventana triangular de 500 puntos con un solapamiento del 50% entre segmentos.

rng default
r = randn(2e4,1);

hx = fir1(30,0.2,rectwin(31));
x = filter(hx,1,r);

hy = ones(1,10);
y = filter(hy,1,r);

win = triang(500);
nov = 250;

cpsd(x,y,win,nov,1024)

Densidades espectrales de potencia cruzada SISO y MIMO

Genere dos sinusoides de dos canales muestreadas a 1 kHz durante 1 segundo. Los canales de la primera señal tienen frecuencias de 200 Hz y 300 Hz. Los canales de la segunda señal tienen frecuencias de 300 Hz y 400 Hz. Ambas señales están integradas en ruido blanco gaussiano de varianza unitaria.

fs = 1e3;
t = (0:1/fs:1-1/fs)';

q = 2*sin(2*pi*[200 300].*t);
q = q+randn(size(q));

r = 2*sin(2*pi*[300 400].*t);
r = r+randn(size(r));

Calcule la densidad espectral de potencia cruzada de las dos señales. Utilice una ventana de Bartlett de 256 muestras para dividir las señales en segmentos y aplique ventanas a los segmentos. Especifique 128 muestras de solapamiento entre segmentos contiguos y 2048 puntos de la DFT. Utilice la funcionalidad incorporada de cpsd para representar el resultado.

cpsd(q,r,bartlett(256),128,2048,fs)

De forma predeterminada, cpsd opera columna a columna en entradas de matriz del mismo tamaño. Cada canal alcanza su pico en las frecuencias de las sinusoides originales.

Repita el cálculo, pero ahora añada 'mimo' a la lista de argumentos.

cpsd(q,r,bartlett(256),128,2048,fs,'mimo')

Cuando se llama con la opción 'mimo', cpsd devuelve un arreglo tridimensional que contiene estimaciones de densidad espectral de potencia cruzada para todas las combinaciones de columnas de entrada. La estimación del segundo canal de q y del primer canal de r muestra un pico realzado en la frecuencia común de 300 Hz.

Fase de espectro cruzado de sinusoides retardadas

Genere dos señales sinusoidales de 100 Hz muestreadas a 1 kHz durante 296 ms. Una de las sinusoides se retrasa 2,5 ms respecto a la otra, lo que equivale a un retraso de π/2. Ambas señales están integradas en ruido blanco gaussiano de varianza 1/4².

Fs = 1000;
t = 0:1/Fs:0.296;

x = cos(2*pi*t*100)+0.25*randn(size(t));
tau = 1/400;
y = cos(2*pi*100*(t-tau))+0.25*randn(size(t));

Calcule y represente la magnitud de la densidad espectral de potencia cruzada. Utilice la configuración predeterminada para cpsd. La magnitud alcanza su pico en la frecuencia en la que existe una coherencia significativa entre las señales.

cpsd(x,y,[],[],[],Fs)

Represente la función de coherencia de magnitud cuadrada y la fase del espectro cruzado. La ordenada en la frecuencia de alta coherencia corresponde al retraso de fase entre las sinusoides.

[Cxy,F] = mscohere(x,y,[],[],[],Fs);

[Pxy,F] = cpsd(x,y,[],[],[],Fs);

subplot(2,1,1)
plot(F,Cxy)
title('Magnitude-Squared Coherence')

subplot(2,1,2)
plot(F,angle(Pxy))

hold on
plot(F,2*pi*100*tau*ones(size(F)),'--')
hold off

xlabel('Hz')
ylabel('\Theta(f)')
title('Cross Spectrum Phase')

Densidad espectral de potencia cruzada de secuencias exponenciales

Genere dos secuencias exponenciales de $$N$$ muestras, $$x_a=a^n$$ y $$x_b=b^n$$, con $$n\ge 0$$. Especifique $$a=0.8$$, $$b=0.9$$ y un $$N$$ pequeño para ver los efectos de tamaño finito.

N = 10;
n = 0:N-1;

a = 0.8;
b = 0.9;

xa = a.^n;
xb = b.^n;

Calcule y represente la densidad espectral de potencia cruzada de las secuencias en el intervalo completo de frecuencias normalizadas, $$[-\pi ,\pi ]$$. Especifique una ventana rectangular de longitud $$N$$ y sin solapamiento entre segmentos.

w = -pi:1/1000:pi;
wind = rectwin(N);
nove = 0;

[pxx,f] = cpsd(xa,xb,wind,nove,w);

El espectro de potencia cruzada de las dos secuencias tiene una expresión analítica para $$N$$ grandes:

$$R(\omega )=\frac{1}{1-a{e}^{-j\omega }}\frac{1}{1-b{e}^{j\omega }}.$$

Convierta esta expresión en una densidad espectral de potencia cruzada dividiéndola por $$2\pi N$$. Compare los resultados. La ondulación en el resultado de cpsd es consecuencia del ventaneo.

nfac = 2*pi*N;

X = 1./(1-a*exp(-1j*w));
Y = 1./(1-b*exp( 1j*w));
R = X.*Y/nfac;

semilogy(f/pi,abs(pxx))
hold on
semilogy(w/pi,abs(R))
hold off
legend('cpsd','Analytic')

Repita el cálculo con $$N=25$$. Las curvas coinciden en seis cifras para $$N$$ tan pequeñas como 100.

N = 25;
n = 0:N-1;
xa = a.^n;
xb = b.^n;

wind = rectwin(N);

[pxx,f] = cpsd(xa,xb,wind,nove,w);
R = X.*Y/(2*pi*N);

semilogy(f/pi,abs(pxx))
hold on
semilogy(w/pi,abs(R))
hold off
legend('cpsd','Analytic')

Reconocimiento del tono de marcación

Utilice la densidad espectral de potencia cruzada para identificar un tono muy dañado.

Las señales sonoras que se generan cuando se marca un número o un símbolo en un teléfono digital son sumas de sinusoides con frecuencias tomadas de dos grupos diferentes. Cada par de tonos contiene una frecuencia del grupo bajo (697 Hz, 770 Hz, 852 Hz o 941 Hz) y una frecuencia del grupo alto (1209 Hz, 1336 Hz o 1477 Hz).

Genere las señales correspondientes a todos los símbolos. Muestree cada tono a 4 kHz durante medio segundo. Prepare una tabla de referencia.

fs = 4e3;
t = 0:1/fs:0.5-1/fs;

nms = ["1";"2";"3";"4";"5";"6";"7";"8";"9";"*";"0";"#"];

ver = [697 770 852 941];
hor = [1209 1336 1477];

v = length(ver);
h = length(hor);

for k = 1:v
    for l = 1:h
        idx = h*(k-1)+l;
        tone = sum(sin(2*pi*[ver(k);hor(l)].*t))';
        tones(:,idx) = tone;
    end
end

Represente el periodograma de Welch de cada señal y anote las frecuencias de los componentes. Utilice una ventana de Hamming de 200 muestras para dividir las señales en segmentos no solapados y aplique ventanas a los segmentos.

[pxx,f] = pwelch(tones,hamming(200),0,[],fs);

for k = 1:v
    for l = 1:h
        idx = h*(k-1)+l;
        ax = subplot(v,h,idx);
        plot(f,10*log10(pxx(:,idx)))
        ylim([-80 0])
        title(nms(idx))
        tx = [ver(k);hor(l)];
        ax.XTick = tx;
        ax.XTickLabel = int2str(tx);
    end
end

La señal que se produce cuando se marca el número 8 se envía a través de un canal ruidoso. La señal recibida está tan dañada que el número no puede identificarse por inspección.

mys = sum(sin(2*pi*[ver(3);hor(2)].*t))'+5*randn(size(t'));

% To hear, type soundsc(mys,fs)

Calcule la densidad espectral de potencia cruzada de la señal dañada y las señales de referencia. Aplique una ventana de Kaiser de 512 muestras con factor de forma β = 5 a las señales. Represente la magnitud de cada espectro.

[pxy,f] = cpsd(mys,tones,kaiser(512,5),100,[],fs);

for k = 1:v
    for l = 1:h
        idx = h*(k-1)+l;
        ax = subplot(v,h,idx);
        plot(f,10*log10(abs(pxy(:,idx))))
        ylim([-80 0])
        title(nms(idx))
        tx = [ver(k);hor(l)];
        ax.XTick = tx;
        ax.XTickLabel = int2str(tx);
    end
end

El dígito de la señal corrupta tiene el espectro con los picos más altos y el valor de RMS más elevado.

[~,loc] = max(rms(abs(pxy)));

digit = nms(loc)

digit = 
"8"

Argumentos de entrada

x, y — Señales de entrada

vectores | matrices

Señales de entrada, especificadas como vectores o matrices.

Ejemplo: cos(pi/4*(0:159))+randn(1,160) especifica una sinusoide integrada en ruido blanco gaussiano.

Tipos de datos: single | double
Soporte de números complejos: Sí

window — Ventana

entero | vector | []

Ventana, especificada como un entero o como un vector fila o vector columna. Utilice window para dividir la señal en segmentos.

• Si window es un entero, cpsd divide x e y en segmentos de longitud window y aplica una ventana de Hamming de esa longitud a cada segmento.

• Si window es un vector, cpsd divide x e y en segmentos de la misma longitud que el vector y aplica ventanas a cada segmento con window.

Si la longitud de x e y no puede dividirse exactamente en un número entero de segmentos con muestras solapadas noverlap, las señales se truncan en consecuencia.

Si se especifica window como vacío, cpsd utiliza una ventana de Hamming de forma que x e y se dividen en ocho segmentos con muestras solapadas noverlap.

Para obtener una lista de las ventanas disponibles, consulte Ventanas.

Ejemplo: hann(N+1) y (1-cos(2*pi*(0:N)'/N))/2 especifican ambas una ventana de Hann de longitud N + 1.

Tipos de datos: single | double

noverlap — Número de muestras solapadas

entero positivo | []

Número de muestras solapadas, especificado como entero positivo.

• Si window es escalar, noverlap debe ser menor que window.

• Si window es un vector, noverlap debe ser menos que la longitud de window.

Si especifica noverlap como vacío, cpsd utiliza un número que genera un 50% de solapamiento entre segmentos. Si no se especifica la longitud del segmento, la función establece noverlap en ⌊N/4,5⌋, donde N es la longitud de las señales de entrada y salida.

Tipos de datos: double | single

nfft — Número de puntos DFT

entero positivo | []

Número de puntos DFT, especificado como entero positivo. Si se especifica nfft como vacío, cpsd establece el parámetro en max(256,2^p), donde p = ⌈log₂ N⌉ para las señales de entrada de longitud N.

Tipos de datos: single | double

fs — Tasa de muestreo

escalar positivo

Tasa de muestreo, especificada como un escalar positivo. La tasa de muestreo es el número de muestras por unidad de tiempo. Si la unidad de tiempo es el segundo, la tasa de muestreo tendrá unidades de Hz.

w — Frecuencias normalizadas

vector

Frecuencias normalizadas, especificadas como vector fila o vector columna con al menos dos elementos. Las frecuencias normalizadas están en rad/muestra.

Ejemplo: w = [pi/4 pi/2]

Tipos de datos: double | single

f — Frecuencias

vector

Frecuencias, especificadas como vector fila o vector columna con al menos dos elementos. Las frecuencias están en ciclos por unidad de tiempo. La unidad de tiempo viene especificada por la tasa de muestreo, fs. Si las unidades de fs son muestras/segundo, las unidades de f son Hz.

Ejemplo: fs = 1000; f = [100 200]

Tipos de datos: double | single

freqrange — Rango de frecuencias para la estimación de la densidad espectral de potencia cruzada

'onesided' | 'twosided' | 'centered'

Rango de frecuencia para la estimación de la densidad espectral de potencia cruzada, especificado como 'onesided', 'twosided' o 'centered'. El valor predeterminado es 'onesided' para señales de valor real y 'twosided' para señales de valor complejo.

• 'onesided': devuelve la estimación unilateral de la densidad espectral de potencia cruzada entre dos señales de entrada de valor real, x e y. Si nfft es par, pxy tiene nfft/2 + 1 filas y se calcula en el intervalo [0,π] rad/muestra. Si nfft es impar, pxy tiene (nfft + 1)/2 filas y el intervalo es [0,π) rad/muestra. Si se especifica fs, los intervalos correspondientes son [0,fs/2] ciclos/unidad de tiempo cuando nfft es par y [0,fs/2) ciclos/unidad de tiempo cuando nfft es impar.

• 'twosided': devuelve la estimación bilateral de la densidad espectral de potencia cruzada de dos señales de entrada de valor real o complejo, x e y. En este caso, pxy tiene nfft filas y se calcula en el intervalo [0,2π) rad/muestra. Si se especifica fs, el intervalo es [0, fs] ciclos/unidad de tiempo.

• 'centered': devuelve la estimación bilateral centrada de la densidad espectral de potencia cruzada de dos señales de entrada de valor real o complejo, x e y. En este caso, pxy tiene nfft filas y se calcula en el intervalo (–π,π] rad/muestra cuando nfft es par y (–π,π) rad/muestra cuando nfft es impar. Si se especifica fs, los intervalos correspondientes son (–fs/2, fs/2] ciclos/unidad de tiempo cuando nfft es par y (–fs/2, fs/2) ciclos/unidad de tiempo cuando nfft es impar.

Si especifica un vector w o f de frecuencias y freqrange como entradas, cpsd ignora freqrange.

Argumentos de salida

pxy — Densidad espectral de potencia cruzada

vector | matriz | arreglo tridimensional

Densidad espectral de potencia cruzada, devuelta como un vector, una matriz o un arreglo tridimensional.

w — Frecuencias normalizadas

vector

Frecuencias normalizadas, devueltas como vector columna de valor real.

• Si pxy es unilateral, w abarca el intervalo [0,π] cuando nfft es par y [0,π) cuando nfft es impar.

• Si pxy es bilateral, w abarca el intervalo [0,2π).

• Si pxy es centrado en DC, w abarca el intervalo (–π,π] cuando nfft es par y (–π,π) cuando nfft es impar.

Tipos de datos: double | single

f — Frecuencias

vector

Frecuencias, devueltas como vector columna de valor real.

Tipos de datos: double | single

Más acerca de

Densidad espectral de potencia cruzada

La densidad espectral de potencia cruzada es la distribución de potencia por unidad de frecuencia y se define como

$$P_{xy}(\omega )={\displaystyle \sum _{m=-\infty }^{\infty }{R}_{xy}(m)e^{-j\omega m}.}$$

La secuencia de correlación cruzada se define como

$$R_{xy}(m)=E\left\{x_{n+m}{y}_n^{\ast }\right\}=E\left\{x_n{y}_{n-m}^{\ast }\right\},$$

donde x_n e y_n son procesos aleatorios conjuntamente estacionarios de media cero, –∞ < n < ∞, $$-\infty <n<\infty$$ y E {· } es el operador de valor esperado.

Algoritmos

cpsd utiliza el método del periodograma modificado y promediado de Welch de la estimación espectral.

Referencias

[1] Oppenheim, Alan V., Ronald W. Schafer, and John R. Buck. Discrete-Time Signal Processing. 2nd Ed. Upper Saddle River, NJ: Prentice Hall, 1999.

[2] Rabiner, Lawrence R., and B. Gold. Theory and Application of Digital Signal Processing. Englewood Cliffs, NJ: Prentice-Hall, 1975, pp. 414–419.

[3] Welch, Peter D. “The Use of the Fast Fourier Transform for the Estimation of Power Spectra: A Method Based on Time Averaging Over Short, Modified Periodograms.” IEEE^® Transactions on Audio and Electroacoustics, Vol. AU-15, June 1967, pp. 70–73.

Capacidades ampliadas

Generación de código C/C++
Genere código C y C++ mediante MATLAB® Coder™.

Arreglos GPU

Acelere código mediante la ejecución en una unidad de procesamiento gráfico (GPU) mediante Parallel Computing Toolbox™.

Esta función es totalmente compatible con los arreglos de GPU. Para obtener más información, consulte Run MATLAB Functions on a GPU (Parallel Computing Toolbox).

Historial de versiones

Introducido antes de R2006a

R2024a: Compatibilidad de la generación de código para entradas de ventana de tamaño variable y precisión simple

La función cpsd admite las entradas de ventana de tamaño variable y precisión simple para la generación de código.

R2023b: Utilizar datos de precisión simple

La función cpsd admite las entradas de precisión simple.

Consulte también

mscohere | pburg | pcov | peig | periodogram | pmcov | pmtm | pmusic | pwelch | pyulear | tfestimate

Temas

• Espectro cruzado y coherencia de magnitud cuadrada