La traducción de esta página aún no se ha actualizado a la versión más reciente. Haga clic aquí para ver la última versión en inglés.

spectrogram

Espectrograma utilizando la transformada de Fourier de tiempo corto

Sintaxis

s = spectrogram(x)

s = spectrogram(x,window)

s = spectrogram(x,window,noverlap)

s = spectrogram(x,window,noverlap,nfft)

[s,w,t] = spectrogram(___)

[s,f,t] = spectrogram(___,fs)

[s,w,t] = spectrogram(x,window,noverlap,w)

[s,f,t] = spectrogram(x,window,noverlap,f,fs)

[___,ps] = spectrogram(___,spectrumtype)

[___] = spectrogram(___,"reassigned")

[___,ps,fc,tc] = spectrogram(___)

[___] = spectrogram(___,freqrange)

[___] = spectrogram(___,Name=Value)

spectrogram(___)

spectrogram(___,freqloc)

Descripción

s = spectrogram(x) devuelve la Transformada de Fourier de tiempo corto (STFT) de la señal de entrada x. Cada columna de s contiene una estimación del contenido de la frecuencia de tiempo corto localizada en el tiempo de x. La magnitud al cuadrado de s se conoce como la representación de tiempo-frecuencia del espectrograma de x [1].

s = spectrogram(x,window) utiliza window para dividir la señal en segmentos y disponerlos en ventanas.

s = spectrogram(x,window,noverlap) utiliza muestras de solapamiento noverlap entre segmentos contiguos.

s = spectrogram(x,window,noverlap,nfft) utiliza puntos de muestreo nfft para calcular la transformada discreta de Fourier.

[s,w,t] = spectrogram(___) devuelve un vector de frecuencias normalizadas, w, y un vector de instantes de tiempo, t, en los que se calcula la STFT. Esta sintaxis puede incluir cualquier combinación de argumentos de entrada de las sintaxis anteriores.

[s,f,t] = spectrogram(___,fs) devuelve un vector de frecuencias cíclicas, f, expresado en términos de la tasa de muestreo fs. fs debe ser la quinta entrada para spectrogram. Para introducir una tasa de muestreo y seguir utilizando los valores predeterminados de los argumentos opcionales anteriores, especifique estos argumentos como vacíos, [].

[s,w,t] = spectrogram(x,window,noverlap,w) devuelve la STFT en las frecuencias normalizadas especificadas en w. w debe tener al menos dos elementos porque, de lo contrario, la función lo interpreta como nfft.

[s,f,t] = spectrogram(x,window,noverlap,f,fs) devuelve la STFT en las frecuencias cíclicas especificadas en f. f debe tener al menos dos elementos porque, de lo contrario, la función lo interpreta como nfft.

[___,ps] = spectrogram(___,spectrumtype) también devuelve una matriz, ps, proporcional al espectrograma de x.

• Si se especifica spectrumtype como "psd", cada columna de ps contiene una estimación de la densidad espectral de potencia (PSD) de un segmento con ventana.

• Si se especifica spectrumtype como "power", cada columna de ps contiene una estimación del espectro de potencia de un segmento con ventana.

[___] = spectrogram(___,"reassigned") reasigna cada estimación de la PSD o del espectro de potencia a la ubicación de su centro de energía. Si su señal contiene componentes temporales o espectrales bien ubicados, esta opción genera un espectrograma más nítido.

[___,ps,fc,tc] = spectrogram(___) también devuelve dos matrices, fc y tc, que contienen la frecuencia y el tiempo del centro de energía de cada estimación de la PSD o del espectro de potencia.

[___] = spectrogram(___,freqrange) devuelve la estimación de la PSD o del espectro de potencia a lo largo del rango de frecuencia especificado por freqrange. Las opciones válidas para freqrange son "onesided", "twosided" y "centered".

[___] = spectrogram(___,Name=Value) especifica opciones adicionales con argumentos de par nombre-valor. Las opciones incluyen el umbral mínimo y la dimensión de tiempo de salida.

spectrogram(___) sin argumentos de salida representa ps en decibelios en la ventana de figura actual.

spectrogram(___,freqloc) especifica el eje en el que se desea representar la frecuencia.

Ejemplos

Valores predeterminados del espectrograma

Genere muestras $$N_x=1024$$ de una señal que consista en una suma de sinusoides. Las frecuencias normalizadas de las sinusoides son $$2\pi /5$$ rad/muestra y $$4\pi /5$$ rad/muestra. La sinusoide de frecuencia superior tiene 10 veces la amplitud de la otra sinusoide.

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

w0 = 2*pi/5;
x = sin(w0*n)+10*sin(2*w0*n);

Calcule la transformada de Fourier de tiempo corto utilizando los valores predeterminados de la función. Represente el espectrograma.

s = spectrogram(x);

spectrogram(x,'yaxis')

Repita el cálculo.

• Divida la señal en secciones de longitud $$nsc=\lfloor N_x/4.5\rfloor$$.

• Disponga en ventanas las secciones utilizando una ventana Hamming.

• Especifique el 50% de solapamiento entre las secciones contiguas.

• Para calcular la FFT, utilice $$\max (256,2^p)$$ puntos, donde $$p=\lceil {\log }_{2}nsc\rceil$$.

Verifique que los dos enfoques dan resultados idénticos.

Nx = length(x);
nsc = floor(Nx/4.5);
nov = floor(nsc/2);
nff = max(256,2^nextpow2(nsc));

t = spectrogram(x,hamming(nsc),nov,nff);

maxerr = max(abs(abs(t(:))-abs(s(:))))

maxerr = 0

Divida la señal en ocho secciones de igual longitud, con un 50% de solapamiento entre las secciones. Especifique la misma longitud de FFT que en el paso anterior. Calcule la transformada de Fourier de tiempo corto y verifique que ofrece el mismo resultado que los dos procedimientos anteriores.

ns = 8;
ov = 0.5;
lsc = floor(Nx/(ns-(ns-1)*ov));

t = spectrogram(x,lsc,floor(ov*lsc),nff);

maxerr = max(abs(abs(t(:))-abs(s(:))))

maxerr = 0

Comparar la función spectrogram y la definición de STFT

Genere una señal compuesta por un chirp cuadrático convexo de valor complejo muestreado a 600 Hz durante 2 segundos. El chirp tiene una frecuencia inicial de 250 Hz y una final de 50 Hz.

fs = 6e2;
ts = 0:1/fs:2;
x = chirp(ts,250,ts(end),50,"quadratic",0,"convex","complex");

Función spectrogram

Utilice la función spectrogram para calcular la STFT de la señal.

• Divida la señal en segmentos, cada uno de $\mathit{M}=49$ muestras.

• Especifique $\mathit{L}=11$ muestras de solapamiento entre segmentos contiguos.

• Descarte el último segmento, el más corto.

• Aplique una ventana de Bartlett a cada segmento.

• Evalúe la transformada discreta de Fourier de cada segmento en ${\mathit{N}}_{\mathrm{DFT}}=1024$ puntos. De forma predeterminada, spectrogram calcula transformadas bilaterales en señales de valor complejo.

M = 49;
L = 11;
g = bartlett(M);
Ndft = 1024;

[s,f,t] = spectrogram(x,g,L,Ndft,fs);

Calcule y muestre el espectrograma, definido como la magnitud al cuadrado de la STFT, con la función waterplot.

waterplot(s,f,t)

Definición de STFT

Calcule la STFT de la señal de ${\mathit{N}}_{\mathit{x}}$ muestras con la definición. Divida la señal en $$\lfloor \frac{N_x-L}{M-L}\rfloor$$ segmentos superpuestos. Aplique una ventana a cada segmento y evalúe la transformada discreta de Fourier en ${\mathit{N}}_{\mathrm{DFT}}$ puntos.

[segs,~] = buffer(1:length(x),M,L,"nodelay");

X = fft(x(segs).*g,Ndft);

Calcule los intervalos de tiempo y frecuencia de la STFT.

• Para encontrar los valores de tiempo, divida el vector tiempo en segmentos superpuestos. Los valores de tiempo son los puntos medios de los segmentos, y cada segmento se trata como un intervalo abierto en el extremo inferior.

• Para encontrar los valores de frecuencia, especifique un intervalo de Nyquist cerrado en la frecuencia cero y abierto en el extremo inferior.

tbuf = ts(segs);
tint = mean(tbuf(2:end,:));

fint = 0:fs/Ndft:fs-fs/Ndft;

Compare la salida de spectrogram con la definición. Muestre el espectrograma.

maxdiff = max(max(abs(s-X)))

maxdiff = 0

waterplot(X,fint,tint)

function waterplot(s,f,t)
% Waterfall plot of spectrogram
    waterfall(f,t,abs(s)'.^2)
    set(gca,XDir="reverse",View=[30 50])
    xlabel("Frequency (Hz)")
    ylabel("Time (s)")
end

Comparar las funciones spectrogram y stft

Genere una señal compuesta por un chirp muestreado a 1,4 kHz durante 2 segundos. La frecuencia del chirp disminuye linealmente de 600 Hz a 100 Hz durante el tiempo de medición.

fs = 1400;
x = chirp(0:1/fs:2,600,2,100);

Valores predeterminados de stft

Calcule la STFT de la señal con las funciones spectrogram y stft. Utilice los valores predeterminados de la función stft:

• Divida la señal en segmentos de 128 muestras y aplique una ventana periódica de Hann a cada segmento.

• Especifique 96 muestras de solapamiento entre segmentos contiguos. Esta longitud equivale al 75% de la longitud de la ventana.

• Especifique 128 puntos de la DFT y centre la STFT en la frecuencia cero, con la frecuencia expresada en hercios.

Compruebe que los dos resultados son iguales.

M = 128;
g = hann(M,"periodic");
L = 0.75*M;
Ndft = 128;

[sp,fp,tp] = spectrogram(x,g,L,Ndft,fs,"centered");

[s,f,t] = stft(x,fs);

dff = max(max(abs(sp-s)))

dff = 0

Represente las dos salidas con la función mesh.

nexttile
mesh(tp,fp,abs(sp).^2)
title("spectrogram")
view(2), axis tight

nexttile
mesh(t,f,abs(s).^2)
title("stft")
view(2), axis tight

Valores predeterminados de spectrogram

Repita el cálculo utilizando los valores predeterminados de la función spectrogram:

• Divida la señal en segmentos de longitud $\mathit{M}=\lfloor {\mathit{N}}_{\mathit{x}}/4.5\rfloor$, donde ${\mathit{N}}_{\mathit{x}}$ es la longitud de la señal. Aplique una ventana de Hamming a cada segmento.

• Especifique el 50% de solapamiento entre segmentos.

• Para calcular la FFT, utilice $\max \left(256,2^{\lceil {\log }_2\mathit{M}\rceil }\right)$ puntos. Calcule el espectrograma solo para las frecuencias normalizadas positivas.

M = floor(length(x)/4.5);
g = hamming(M);
L = floor(M/2);
Ndft = max(256,2^nextpow2(M));

[sx,fx,tx] = spectrogram(x);

[st,ft,tt] = stft(x,Window=g,OverlapLength=L, ...
    FFTLength=Ndft,FrequencyRange="onesided");

dff = max(max(sx-st))

dff = 0

Represente las dos salidas con la función waterplot. Divida el eje de la frecuencia por $\pi$ en ambos casos. Para la salida de stft, divida los números de muestra por la tasa de muestreo efectiva, $2\pi$.

figure
nexttile
waterplot(sx,fx/pi,tx)
title("spectrogram")

nexttile
waterplot(st,ft/pi,tt/(2*pi))
title("stft")

function waterplot(s,f,t)
% Waterfall plot of spectrogram
    waterfall(f,t,abs(s)'.^2)
    set(gca,XDir="reverse",View=[30 50])
    xlabel("Frequency/\pi")
    ylabel("Samples")
end

Espectrograma y frecuencia instantánea

Utilice la función spectrogram para medir y realizar un seguimiento de la frecuencia instantánea de una señal.

Genere un chirp cuadrático muestreado a 1 kHz durante dos segundos. Especifique el chirp de forma que su frecuencia sea inicialmente 100 Hz y aumente a 200 Hz al pasar un segundo.

fs = 1000;
t = 0:1/fs:2-1/fs;
y = chirp(t,100,1,200,'quadratic');

Calcule el espectro del chirp utilizando la transformada de Fourier de tiempo corto implementada en la función spectrogram. Divida la señal en secciones de longitud 100, dispuestas en ventanas con una ventana Hamming. Especifique 80 muestras de solapamiento entre secciones contiguas y evalúe el espectro en frecuencias $$\lfloor 100/2+1\rfloor =51$$.

spectrogram(y,100,80,100,fs,'yaxis')

Realice un seguimiento de la frecuencia del chirp buscando la cresta de tiempo-frecuencia con más energía en los puntos de tiempo $$\lfloor (2000-80)/(100-80)\rfloor =96$$. Superponga la frecuencia instantánea en la gráfica del espectrograma.

[~,f,t,p] = spectrogram(y,100,80,100,fs);

[fridge,~,lr] = tfridge(p,f);

hold on
plot3(t,fridge,abs(p(lr)),'LineWidth',4)
hold off

Espectrograma de señal compleja

Genere 512 muestras de un chirp con contenido de frecuencia que varíe sinusoidalmente.

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

x = exp(1j*pi*sin(8*n/N)*32);

Calcule la transformada de Fourier de tiempo corto bilateral centrada del chirp. Divida la señal en segmentos de 32 muestras con un solapamiento de 16 muestras. Especifique 64 puntos DFT. Represente el espectrograma.

[scalar,fs,ts] = spectrogram(x,32,16,64,'centered');

spectrogram(x,32,16,64,'centered','yaxis')

Obtenga el mismo resultado calculando el espectrograma en 64 frecuencias equiespaciadas a lo largo del intervalo $(-\pi ,\pi ]$. La opción 'centered' no es necesaria.

fintv = -pi+pi/32:pi/32:pi;

[vector,fv,tv] = spectrogram(x,32,16,fintv);

spectrogram(x,32,16,fintv,'yaxis')

Comparar las funciones spectrogram y pspectrum

Genere una señal que conste de un oscilador controlado por tensión y tres átomos gaussianos. La señal se muestrea a ${\mathit{f}}_{\mathit{s}}=2$ kHz durante 2 segundos.

fs = 2000;
tx = 0:1/fs:2;
gaussFun = @(A,x,mu,f) exp(-(x-mu).^2/(2*0.03^2)).*sin(2*pi*f.*x)*A';
s = gaussFun([1 1 1],tx',[0.1 0.65 1],[2 6 2]*100)*1.5;
x = vco(chirp(tx+.1,0,tx(end),3).*exp(-2*(tx-1).^2),[0.1 0.4]*fs,fs);
x = s+x';

Transformadas de Fourier de tiempo corto

Utilice la función pspectrum para calcular la STFT.

• Divida la señal de ${\mathit{N}}_{\mathit{x}}$ muestras en segmentos con longitud $\mathit{M}=80$ muestras, correspondientes a una resolución temporal de $80/2000=40$ milisegundos.

• Especifique $\mathit{L}=16$ muestras o el 20% de solapamiento entre los segmentos contiguos.

• Añada una ventana de Kaiser a cada segmento y especifique un manchado de $\ell =0.7$.

M = 80;
L = 16;
lk = 0.7;

[S,F,T] = pspectrum(x,fs,"spectrogram", ...
    TimeResolution=M/fs,OverlapPercent=L/M*100, ...
    Leakage=lk);

Compárela con el resultado obtenido con la función spectrogram.

• Especifique la longitud de la ventana y el solapamiento directamente en las muestras.

• pspectrum siempre usa una ventana de Kaiser como $\mathit{g}\left(\mathit{n}\right)$. El manchado $\ell$ y el factor de forma $\beta$ de la ventana están relacionados por $\beta =40\times \left(1-\ell \right)$.

• pspectrum siempre usa ${\mathit{N}}_{\mathrm{DFT}}=1024$ puntos cuando calcula la transformada discreta de Fourier. Puede especificar este número si desea calcular la transformada en un rango de frecuencias bilateral o centrado. Sin embargo, en las transformadas unilaterales, que son las predeterminadas para las señales reales, spectrogram utiliza $1024/2+1=513$ puntos. Como alternativa, puede especificar el vector de frecuencias en el que desea calcular la transformada, como en este ejemplo.

• Si una señal no puede dividirse exactamente en $\mathit{k}=\lfloor \frac{{\mathit{N}}_{\mathit{x}}-\mathit{L}}{\mathit{M}-\mathit{L}}\rfloor$ segmentos, spectrogram la trunca, mientras que pspectrum la rellena con ceros para crear un segmento extra. Para que las salidas sean equivalentes, elimine el segmento final y el elemento final del vector de tiempo.

• spectrogram devuelve la STFT, cuya magnitud elevada al cuadrado es el espectrograma. pspectrum devuelve el espectro de potencia segmento a segmento, que ya está elevado al cuadrado pero se divide por un factor de ${\sum }_{\mathit{n}}\mathit{g}\left(\mathit{n}\right)$ antes de elevarlo al cuadrado.

• En las transformadas unilaterales, pspectrum añade un factor extra de 2 al espectrograma.

g = kaiser(M,40*(1-lk));

k = (length(x)-L)/(M-L);
if k~=floor(k)
    S = S(:,1:floor(k));
    T = T(1:floor(k));
end

[s,f,t] = spectrogram(x/sum(g)*sqrt(2),g,L,F,fs);

Para mostrar los espectrogramas calculados por las dos funciones, utilice la función waterplot.

subplot(2,1,1)
waterplot(sqrt(S),F,T)
title("pspectrum")

subplot(2,1,2)
waterplot(s,f,t)
title("spectrogram")

maxd = max(max(abs(abs(s).^2-S)))

maxd = 2.4419e-08

Espectros de potencia y gráficas de conveniencia

La función spectrogram tiene un cuarto argumento que se corresponde con el espectro de potencia segmento a segmento o con la densidad espectral de potencia. Al igual que la salida de pspectrum, el argumento ps ya está elevado al cuadrado e incluye el factor de normalización ${\sum }_{\mathit{n}}\mathit{g}\left(\mathit{n}\right)$. En los espectrogramas unilaterales de señales reales, hay que incluir el factor adicional de 2. Establezca el argumento de escala de la función en "power".

[~,~,~,ps] = spectrogram(x*sqrt(2),g,L,F,fs,"power");

max(abs(S(:)-ps(:)))

ans = 2.4419e-08

Cuando se llaman sin argumentos de salida, tanto pspectrum como spectrogram representan el espectrograma de la señal en decibelios. Incluya el factor de 2 en los espectrogramas unilaterales. Establezca que los mapas de colores sean los mismos en ambas gráficas. Fije los límites del eje x en los mismos valores para hacer visible el segmento extra al final de la gráfica de pspectrum. En la gráfica de spectrogram, muestre la frecuencia en el eje y.

subplot(2,1,1)
pspectrum(x,fs,"spectrogram", ...
    TimeResolution=M/fs,OverlapPercent=L/M*100, ...
    Leakage=lk)
title("pspectrum")
cc = clim;
xl = xlim;

subplot(2,1,2)
spectrogram(x*sqrt(2),g,L,F,fs,"power","yaxis")
title("spectrogram")
clim(cc)
xlim(xl)

function waterplot(s,f,t)
% Waterfall plot of spectrogram
    waterfall(f,t,abs(s)'.^2)
    set(gca,XDir="reverse",View=[30 50])
    xlabel("Frequency (Hz)")
    ylabel("Time (s)")
end

Espectrograma reasignado de chirp cuadrático

Genere una señal de chirp muestreada durante 2 segundos a 1 Hz. Especifique el chirp de forma que su frecuencia sea inicialmente 100 Hz y aumente a 200 Hz al pasar 1 segundo.

fs = 1000;
t = 0:1/fs:2;
y = chirp(t,100,1,200,"quadratic");

Calcule el espectrograma reasignado de la señal.

• Divida la señal en secciones de longitud 128, dispuestas en ventanas con una ventana de Kaiser con el parámetro de forma $$\beta =18$$.

• Especifique 120 muestras de solapamiento entre secciones contiguas.

• Evalúe el espectro en frecuencias $$\lfloor 128/2\rfloor =65$$ y bins de tiempo $$\lfloor (length(x)-120)/(128-120)\rfloor =235$$.

Utilice la función spectrogram sin argumentos de salida para representar el espectrograma reasignado. Muestre la frecuencia en el eje y y el tiempo en el eje x.

spectrogram(y,kaiser(128,18),120,128,fs, ...
    "reassigned","yaxis")

Vuelva a generar la gráfica usando la función imagesc. Especifique la dirección del eje y para que los valores de frecuencia aumenten de abajo hacia arriba. Añada eps al espectrograma reasignado para evitar posibles infinitos negativos al convertir a decibelios.

[~,fr,tr,pxx] = spectrogram(y,kaiser(128,18),120,128,fs, ...
    "reassigned");

imagesc(tr,fr,pow2db(pxx+eps))
axis xy
xlabel("Time (s)")
ylabel("Frequency (Hz)")
colorbar

Espectrograma con umbral

Genere una señal de chirp muestreada durante 2 segundos a 1 Hz. Especifique el chirp de forma que su frecuencia sea inicialmente 100 Hz y aumente a 200 Hz al pasar 1 segundo.

Fs = 1000;
t = 0:1/Fs:2;
y = chirp(t,100,1,200,'quadratic');

Calcule la densidad espectral de potencia (PSD) dependiente del tiempo de la señal.

• Divida la señal en secciones de longitud 128, dispuestas en ventanas con una ventana de Kaiser con el parámetro de forma $$\beta =18$$.

• Especifique 120 muestras de solapamiento entre secciones contiguas.

• Evalúe el espectro en frecuencias $$\lfloor 128/2\rfloor =65$$ y bins de tiempo $$\lfloor (length(x)-120)/(128-120)\rfloor =235$$.

Obtenga como salida la frecuencia y el tiempo del centro de gravedad de cada estimación de PSD. Establezca en cero los elementos de la PSD inferiores a $$-30$$ dB.

[~,~,~,pxx,fc,tc] = spectrogram(y,kaiser(128,18),120,128,Fs, ...
    'MinThreshold',-30);

Represente los elementos distintos de cero como funciones de las frecuencias y tiempos del centro de gravedad.

plot(tc(pxx>0),fc(pxx>0),'.')

Cálculo de espectrogramas centrados y unilaterales

Genere una señal compuesta por un chirp de valor real muestreado a 2 kHz durante 2 segundos.

fs = 2000;
tx = 0:1/fs:2;
x = vco(-chirp(tx,0,tx(end),2).*exp(-3*(tx-1).^2), ...
    [0.1 0.4]*fs,fs).*hann(length(tx))';

Espectrograma bilateral

Calcule y represente la STFT bilateral de la señal.

• Divida la señal en segmentos, cada uno de $\mathit{M}=73$ muestras.

• Especifique $\mathit{L}=24$ muestras de solapamiento entre segmentos contiguos.

• Descarte el último segmento, el más corto.

• Aplique una ventana flat-top a cada segmento.

• Evalúe la transformada discreta de Fourier de cada segmento en ${\mathit{N}}_{\mathrm{DFT}}=895$ puntos, teniendo en cuenta que es un número impar.

M = 73;
L = 24;
g = flattopwin(M);
Ndft = 895;
neven = ~mod(Ndft,2);

[stwo,f,t] = spectrogram(x,g,L,Ndft,fs,"twosided");

Utilice la función spectrogram sin argumentos de salida para representar el espectrograma bilateral.

spectrogram(x,g,L,Ndft,fs,"twosided","power","yaxis");

Calcule el espectrograma bilateral con la definición. Divida la señal en segmentos de $\mathit{M}$ muestras con $\mathit{L}$ muestras de solapamiento entre segmentos contiguos. Aplique una ventana a cada segmento y calcule la transformada discreta de Fourier en ${\mathit{N}}_{\mathrm{DFT}}$ puntos.

[segs,~] = buffer(1:length(x),M,L,"nodelay");

Xtwo = fft(x(segs).*g,Ndft);

Calcule los intervalos de tiempo y frecuencia.

• Para encontrar los valores de tiempo, divida el vector tiempo en segmentos superpuestos. Los valores de tiempo son los puntos medios de los segmentos, y cada segmento se trata como un intervalo abierto en el extremo inferior.

• Para encontrar los valores de frecuencia, especifique un intervalo de Nyquist cerrado en la frecuencia cero y abierto en el extremo superior.

tbuf = tx(segs);
ttwo = mean(tbuf(2:end,:));

ftwo = 0:fs/Ndft:fs*(1-1/Ndft);

Compare las salidas del spectrogram con las definiciones. Visualice los espectrogramas con la función waterplot.

diffs = [max(max(abs(stwo-Xtwo))) max(abs(f-ftwo')) max(abs(t-ttwo))]

diffs = 1×3
10-12 ×

         0    0.2274    0.0002

figure
nexttile
waterplot(Xtwo,ftwo,ttwo)
title("Two-Sided, Definition")

nexttile
waterplot(stwo,f,t)
title("Two-Sided, spectrogram Function")

Espectrograma centrado

Calcule el espectrograma centrado de la señal.

• Utilice los mismos valores de tiempo que utilizó en la STFT bilateral.

• Desplace el componente de frecuencia cero de la STFT hacia el centro del espectro con la función fftshift.

• En los valores impares de ${\mathit{N}}_{\mathrm{DFT}}$, el intervalo de frecuencia está abierto en ambos extremos. En los valores pares de ${\mathit{N}}_{\mathrm{DFT}}$, el intervalo de frecuencia está abierto en el extremo inferior y cerrado en el extremo superior.

Compare las salidas y muestre los espectrogramas.

tcen = ttwo;

if ~neven
    Xcen = fftshift(Xtwo,1);
    fcen = -fs/2*(1-1/Ndft):fs/Ndft:fs/2;
else
    Xcen = fftshift(circshift(Xtwo,-1),1);
    fcen = (-fs/2*(1-1/Ndft):fs/Ndft:fs/2)+fs/Ndft/2;
end

[scen,f,t] = spectrogram(x,g,L,Ndft,fs,"centered");

diffs = [max(max(abs(scen-Xcen))) max(abs(f-fcen')) max(abs(t-tcen))]

diffs = 1×3
10-12 ×

         0    0.2274    0.0002

figure
nexttile
waterplot(Xcen,fcen,tcen)
title("Centered, Definition")

nexttile
waterplot(scen,f,t)
title("Centered, spectrogram Function")

Espectrograma unilateral

Calcule el espectrograma unilateral de la señal.

• Utilice los mismos valores de tiempo que utilizó en la STFT bilateral.

• En los valores impares de ${\mathit{N}}_{\mathrm{DFT}}$, la STFT unilateral consiste en las primeras $\left({\mathit{N}}_{\mathrm{DFT}}+1\right)/2$ filas de la STFT bilateral. En los valores pares de ${\mathit{N}}_{\mathrm{DFT}}$, la STFT unilateral consiste en las primeras ${\mathit{N}}_{\mathrm{DFT}}/2+1$ filas de la STFT bilateral.

• En los valores impares de ${\mathit{N}}_{\mathrm{DFT}}$, el intervalo de frecuencia se cierra en la frecuencia cero y se abre en la frecuencia de Nyquist. En los valores pares de ${\mathit{N}}_{\mathrm{DFT}}$, el intervalo de frecuencia se cierra en ambos extremos.

Compare las salidas y muestre los espectrogramas. En señales de valor real, el argumento "onesided" es opcional.

tone = ttwo;

if ~neven
    Xone = Xtwo(1:(Ndft+1)/2,:);
else
    Xone = Xtwo(1:Ndft/2+1,:);
end

fone = 0:fs/Ndft:fs/2;

[sone,f,t] = spectrogram(x,g,L,Ndft,fs);

diffs = [max(max(abs(sone-Xone))) max(abs(f-fone')) max(abs(t-tone))]

diffs = 1×3
10-12 ×

         0    0.1137    0.0002

figure
nexttile
waterplot(Xone,fone,tone)
title("One-Sided, Definition")

nexttile
waterplot(sone,f,t)
title("One-Sided, spectrogram Function")

function waterplot(s,f,t)
% Waterfall plot of spectrogram
waterfall(f,t,abs(s)'.^2)
set(gca,XDir="reverse",View=[30 50])
xlabel("Frequency (Hz)")
ylabel("Time (s)")
end

Calcular las PSD de los segmentos y los espectros de potencia

La función spectrogram tiene como cuarto argumento de salida una matriz que contiene la densidad espectral de potencia (PSD) o el espectro de potencia de cada segmento. El espectro de potencia es igual a la PSD multiplicada por el ancho de banda de ruido equivalente (ENBW) de la ventana.

Genere una señal compuesta por un chirp logarítmico muestreado a 1 kHz durante 1 segundo. El chirp tiene una frecuencia inicial de 400 Hz que disminuye a 10 Hz al final de la medición.

fs = 1000;
tt = 0:1/fs:1-1/fs;
y = chirp(tt,400,tt(end),10,"logarithmic");

PSD de segmentos y espectros de potencia con tasa de muestreo

Divida la señal en segmentos de 102 muestras y aplique una ventana de Hann a cada uno. Especifique 12 muestras de solapamiento entre segmentos contiguos y 1024 puntos de la DFT.

M = 102;
g = hann(M);
L = 12;
Ndft = 1024;

Calcule el espectrograma de la señal con el tipo de espectro de la PSD predeterminado. Genere la STFT y el arreglo de densidades espectrales de potencia de los segmentos.

[s,f,t,p] = spectrogram(y,g,L,Ndft,fs);

Repita el cálculo con el tipo de espectro especificado como "power". Obtenga como salida la STFT y el arreglo de espectros de potencia de los segmentos.

[r,~,~,q] = spectrogram(y,g,L,Ndft,fs,"power");

Compruebe que el espectrograma es el mismo en ambos casos. Represente el espectrograma con una escala logarítmica para la frecuencia.

max(max(abs(s).^2-abs(r).^2))

ans = 0

waterfall(f,t,abs(s)'.^2)
set(gca,XScale="log",...
    XDir="reverse",View=[30 50])

Compruebe que los espectros de potencia son iguales a las densidades espectrales de potencia multiplicadas por el ENBW de la ventana.

max(max(abs(q-p*enbw(g,fs))))

ans = 1.1102e-16

Compruebe que la matriz de espectros de potencia del segmento es proporcional al espectrograma. El factor de proporcionalidad es el cuadrado de la suma de los elementos de la ventana.

max(max(abs(s).^2-q*sum(g)^2))

ans = 3.4694e-18

PSD de segmentos y espectros de potencia con frecuencias normalizadas

Repita el cálculo, pero ahora trabaje en frecuencias normalizadas. Los resultados son los mismos cuando se especifica la tasa de muestreo como $2\pi$.

[~,~,~,pn] = spectrogram(y,g,L,Ndft);
[~,~,~,qn] = spectrogram(y,g,L,Ndft,"power");

max(max(abs(qn-pn*enbw(g,2*pi))))

ans = 1.1102e-16

Seguimiento de chirps en una señal de audio

Cargue una señal de audio que contenga dos chirps decrecientes y un sonido de salpicadura de banda ancha. Calcule la transformada de Fourier de tiempo corto. Divida la forma de onda en segmentos de 400 muestras con 300 muestras de solapamiento. Represente el espectrograma.

load splat

% To hear, type soundsc(y,Fs)

sg = 400;
ov = 300;

spectrogram(y,sg,ov,[],Fs,"yaxis")
colormap bone

Utilice la función spectrogram para emitir la densidad espectral de potencia (PSD) de la señal.

[s,f,t,p] = spectrogram(y,sg,ov,[],Fs);

Realice un seguimiento de los dos chirps utilizando la función medfreq. Para encontrar el chirp de baja frecuencia más potente, restrinja la búsqueda a las frecuencias por encima de 100 Hz y a los tiempos anteriores al inicio del sonido de banda ancha.

f1 = f > 100;
t1 = t < 0.75;

m1 = medfreq(p(f1,t1),f(f1));

Para encontrar el chirp de baja frecuencia más débil, restrinja la búsqueda a las frecuencias por encima de 2500 Hz y a tiempos entre 0,3 segundos y 0,65 segundos.

f2 = f > 2500;
t2 = t > 0.3 & t < 0.65;

m2 = medfreq(p(f2,t2),f(f2));

Superponga el resultado en el espectrograma. Divida los valores de frecuencia entre 1000 para expresarlos en kHz.

hold on
plot(t(t1),m1/1000,LineWidth=4)
plot(t(t2),m2/1000,LineWidth=4)
hold off

Visualización de espectrogramas 3D

Genere dos segundos de una señal muestreada a 10 kHz. Especifique la frecuencia instantánea de la señal como una función triangular de tiempo.

fs = 10e3;
t = 0:1/fs:2;
x1 = vco(sawtooth(2*pi*t,0.5),[0.1 0.4]*fs,fs);

Calcule y represente el espectrograma de la señal. Utilice una ventana de Kaiser de longitud 256 y parámetro de forma $$\beta =5$$. Especifique 220 muestras de solapamiento sección a sección y 512 puntos DFT. Represente la frecuencia en el eje y. Utilice el mapa de colores y la vista predeterminados.

spectrogram(x1,kaiser(256,5),220,512,fs,'yaxis')

Cambie la vista para mostrar el espectrograma como gráfica de cascada. Establezca el mapa de colores en bone.

view(-45,65)
colormap bone

Argumentos de entrada

x — Señal de entrada

vector

Señal de entrada, especificada como un vector fila o vector columna.

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, spectrogram divide x en segmentos de longitud window y dispone en ventanas cada segmento con una ventana Hamming de esa longitud.

• Si window es un vector, spectrogram divide x en segmentos de la misma longitud que el vector y dispone en ventanas cada segmento utilizando window.

Si la longitud de x no puede dividirse exactamente en un número entero de segmentos con muestras solapadas noverlap, x se trunca en consecuencia.

Si especifica window como vacío, spectrogram utiliza una ventana Hamming de forma que x se divide 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.

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, spectrogram 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_x/4.5⌋, donde N_x es la longitud de la señal de entrada y los símbolos ⌊⌋ determinan la función de suelo.

nfft — Número de puntos DFT

escalar entero positivo | []

Número de puntos DFT, especificado como escalar entero positivo. Si especifica nfft como vacío, spectrogram establece el parámetro en max(256,2^p), donde p = ⌈log₂ N_w⌉, los símbolos ⌈⌉ determinan la función de techo, y

• N_w = window si window es un escalar.

• N_w = length(window) si window es un vector.

w — Frecuencias normalizadas

vector

Frecuencias normalizadas, especificadas como vector. w debe tener al menos dos elementos porque, de lo contrario, la función lo interpreta como nfft. Las frecuencias normalizadas están en rad/muestra.

Ejemplo: pi./[2 4]

f — Frecuencias cíclicas

vector

Frecuencias cíclicas, especificadas como vector. f debe tener al menos dos elementos porque, de lo contrario, la función lo interpreta como nfft. Las unidades de f vienen especificadas por la tasa de muestreo, fs.

fs — Tasa de muestreo

1 Hz (predeterminado) | 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 estará en Hz.

freqrange — Rango de frecuencia para estimación de PSD

"onesided" | "twosided" | "centered"

Rango de frecuencia para la estimación de la PSD, especificado como "onesided", "twosided" o "centered". En señales de valor real, el valor predeterminado es "onesided". En señales de valor complejo, el valor predeterminado es "twosided", y especificar "onesided" da como resultado un error.

• "onesided" devuelve el espectrograma unilateral de una señal de entrada real. Si nfft es par, ps tiene nfft/2 + 1 filas y se calcula en el intervalo [0, π] rad/muestra. Si nfft es impar, ps tiene (nfft + 1)/2 filas y el intervalo es [0, π) rad/muestra. Si especifica fs, los intervalos son respectivamente [0, fs/2] ciclos/unidad de tiempo y [0, fs/2) ciclos/unidad de tiempo.

• "twosided" devuelve el espectrograma bilateral de una señal de valor real o complejo. ps tiene nfft filas y se calcula en el intervalo [0, 2π) rad/muestra. Si especifica fs, el intervalo es [0, fs] ciclos/unidad de tiempo.

• "centered" devuelve el espectrograma bilateral centrado de una señal de valor real o complejo. ps tiene nfft filas. Si nfft es par, ps se calcula en el intervalo (–π, π] rad/muestra. Si nfft es impar, ps se calcula en (–π, π) rad/muestra. Si especifica fs, los intervalos son respectivamente (–fs/2, fs/2] ciclos/unidad de tiempo y (–fs/2, fs/2) ciclos/unidad de tiempo.

spectrumtype — Escalado del espectro de potencia

"psd" (predeterminado) | "power"

Escalado del espectro de potencia, especificado como "psd" o "power".

• Omitir spectrumtype o especificar "psd" devuelve la densidad espectral de potencia.

• Especificando "power", se escala cada estimación de PSD por el ancho de banda de ruido equivalente de la ventana. El resultado es una estimación de la potencia en cada frecuencia. Si la opción "reassigned" está activada, la función integra la PSD en la anchura de cada bin de frecuencia antes de reasignar.

freqloc — Eje de visualización de la frecuencia

"xaxis" (predeterminado) | "yaxis"

Eje de visualización de la frecuencia, especificado como "xaxis" o "yaxis".

• "xaxis" muestra la frecuencia en el eje x y el tiempo en el eje y.

• "yaxis" muestra la frecuencia en el eje y y el tiempo en el eje x.

Este argumento se ignora si se llama a spectrogram con argumentos de salida.

Argumentos de par nombre-valor

Especifique pares de argumentos opcionales como Name1=Value1,...,NameN=ValueN, donde Name es el nombre del argumento y Value es el valor correspondiente. Los argumentos de nombre-valor deben aparecer después de otros argumentos. Sin embargo, el orden de los pares no importa.

Ejemplo: spectrogram(x,100,OutputTimeDimension="downrows") divide x en segmentos de longitud 100 y dispone en ventanas cada segmento con una ventana Hamming de esa longitud. La salida del espectrograma tiene una dimensión de tiempo en las filas.

En las versiones anteriores a la R2021a, utilice comas para separar cada nombre y valor, y encierre Name entre comillas.

Ejemplo: spectrogram(x,100,'OutputTimeDimension','downrows') divide x en segmentos de longitud 100 y dispone en ventanas cada segmento con una ventana Hamming de esa longitud. La salida del espectrograma tiene una dimensión de tiempo en las filas.

MinThreshold — Umbral

-Inf (predeterminado) | escalar real

Umbral, especificado como un escalar real expresado en decibelios. spectrogram establece en cero esos elementos de s de forma que 10 log₁₀(s) ≤ thresh.

OutputTimeDimension — Dimensión de tiempo de salida

"acrosscolumns" (predeterminado) | "downrows"

Dimensión de tiempo de salida, especificada como "acrosscolumns" o "downrows". Establezca este valor en "downrows" si desea la dimensión de tiempo de s, ps, fc y tc en las filas y la dimensión de frecuencia en las columnas. Establezca este valor en "acrosscolumns" si desea la dimensión de tiempo de s, ps, fc y tc en las columnas y la dimensión de frecuencia en las filas. Esta entrada se ignora si se llama a la función sin argumentos de salida.

Argumentos de salida

s — Transformada de Fourier de tiempo corto

matriz

Transformada de Fourier de tiempo corto, devuelta como matriz. El tiempo aumenta en las columnas de s y la frecuencia se reduce en las filas, empezando desde cero.

• Si x es una señal de longitud N_x, s tiene columnas k, donde

  • k = ⌊(N_x – noverlap)/(window – noverlap)⌋ si window es un escalar.

  • k = ⌊(N_x – noverlap)/(length(window) – noverlap)⌋ si window es un vector.

• Si x es real y nfft es par, s tiene (nfft/2 + 1) filas.

• Si x es real y nfft es impar, s tiene (nfft + 1)/2 filas.

• Si x es de valor complejo, s tiene nfft filas.

Nota

Cuando freqrange se establece en "onesided", spectrogram emite los valores s en el rango positivo de Nyquist y no conserva la potencia total.

s no se ve afectada por la opción "reassigned".

w — Frecuencias normalizadas

vector

Frecuencias normalizadas, devueltas como vector. w tiene una longitud igual al número de filas de s.

t — Instantes de tiempo

vector

Instantes de tiempo, devueltos como vector. Los valores de tiempo en t corresponden al punto de en medio de cada segmento.

f — Frecuencias cíclicas

vector

Frecuencias cíclicas, devueltas como vector. f tiene una longitud igual al número de filas de s.

ps — Densidad espectral de potencia o espectro de potencia

matriz

Densidad espectral de potencia (PSD) o espectro de potencia, devueltos como matriz.

• Si x es real y freqrange se deja sin especificar o se establece en "onesided", ps contiene la estimación del periodograma unilateral modificada de la PSD o del espectro de potencia de cada segmento. La función multiplica la potencia por 2 en todas las frecuencias excepto 0 y la frecuencia de Nyquist para conservar la potencia total.

• Si x es de valor complejo o si freqrange se establece en "twosided" o "centered", ps contiene la estimación del periodograma bilateral modificada de la PSD o del espectro de potencia de cada segmento.

• Si se especifica un vector de frecuencias normalizadas en w o un vector de frecuencias cíclicas en f, ps contiene la estimación del periodograma modificada de la PSD o espectro de potencia de cada segmento evaluado en las frecuencias de entrada.

fc, tc — Frecuencias y tiempos del centro de energía

matrices

Frecuencias y tiempos del centro de energía, devueltos como matrices del mismo tamaño que la transformada de Fourier de tiempo corto. Si no especifica una tasa de muestreo, los elementos de fc se devuelven como frecuencias normalizadas.

Más acerca de

Transformada de Fourier de tiempo corto

La transformada de Fourier de tiempo corto (STFT) se utiliza para analizar cómo cambia el contenido de frecuencia de una señal no estacionaria a lo largo del tiempo. La magnitud al cuadrado de la STFT se conoce como la representación de tiempo-frecuencia del espectrograma de la señal. Para obtener más información sobre el espectrograma y cómo calcularlo mediante las funciones de Signal Processing Toolbox™, consulte Spectrogram Computation with Signal Processing Toolbox.

La STFT de una señal se calcula deslizando una ventana de análisis g(n) de longitud M sobre la señal y calculando la transformada discreta de Fourier (DFT) de cada segmento de los datos con ventana. La ventana salta sobre la señal original en intervalos de R muestras, lo que equivale a L = M – R muestras de solapamiento entre segmentos contiguos. La mayoría de las funciones de ventana se estrechan en los bordes para evitar el anillamiento espectral. La DFT de cada segmento con ventana se añade a una matriz de valor complejo que contiene la magnitud y la fase de cada punto en tiempo y frecuencia. La matriz STFT tiene

$$k=\lfloor \frac{N_x-L}{M-L}\rfloor$$

 columnas, donde N_x es la longitud de la señal x(n) y los símbolos ⌊⌋ denotan la función de suelo. En las transformadas centradas y bilaterales, el número de filas de la matriz es igual a N_DFT, el número de puntos de la DFT, y en las transformadas unilaterales de señales de valor real es un número impar cercano a N_DFT/2.

La columna m de la matriz de la STFT $$X(f)=\left[\begin{array}{ccccc}{X}_1(f)& X_2(f)& X_3(f)& \cdots & X_k(f)\end{array}\right]$$ contiene la DFT de los datos con ventana centrados en el tiempo mR:

$$X_m(f)={\displaystyle \sum _{n=-\infty }^{\infty }x(n)g(n-mR)e^{-j2\pi fn}}.$$

• La transformada de Fourier de tiempo corto es invertible. El proceso de inversión solapa y suma los segmentos con ventana para compensar la atenuación de la señal en los bordes de la ventana. Para obtener más información, consulte Inverse Short-Time Fourier Transform.

• La función istft invierte la STFT de una señal.

• En una serie de circunstancias concretas es posible lograr la "reconstrucción perfecta" de una señal. Para obtener más información, consulte Perfect Reconstruction.

• La función stftmag2sig devuelve una estimación de una señal reconstruida a partir de la magnitud de su STFT.

Sugerencias

Si una transformada de Fourier de tiempo corto tiene ceros, su conversión a decibelios da como resultado infinitos negativos que no pueden representarse. Para evitar esta posible dificultad, spectrogram añade eps a la transformada de Fourier de tiempo corto cuando la llama sin argumentos de salida.

Referencias

[1] Boashash, Boualem, ed. Time Frequency Signal Analysis and Processing: A Comprehensive Reference. Second edition. EURASIP and Academic Press Series in Signal and Image Processing. Amsterdam and Boston: Academic Press, 2016.

[2] Chassande-Motin, Éric, François Auger, and Patrick Flandrin. "Reassignment." In Time-Frequency Analysis: Concepts and Methods. Edited by Franz Hlawatsch and François Auger. London: ISTE/John Wiley and Sons, 2008.

[3] Fulop, Sean A., and Kelly Fitz. "Algorithms for computing the time-corrected instantaneous frequency (reassigned) spectrogram, with applications." Journal of the Acoustical Society of America. Vol. 119, January 2006, pp. 360–371.

[4] Oppenheim, Alan V., and Ronald W. Schafer, with John R. Buck. Discrete-Time Signal Processing. Second edition. Upper Saddle River, NJ: Prentice Hall, 1999.

[5] Rabiner, Lawrence R., and Ronald W. Schafer. Digital Processing of Speech Signals. Englewood Cliffs, NJ: Prentice-Hall, 1978.

Capacidades ampliadas

Arreglos altos

Realice cálculos con arreglos que tienen más filas de las que caben en la memoria.

Notas y limitaciones de uso:

• La entrada debe ser un vector columna alto.

• El argumento window debe especificarse siempre.

• OutputTimeDimension debe especificarse siempre y establecerse en "downrows".

• La opción reassigned no es compatible.

• Las sintaxis sin argumentos de salida no son compatibles.

Para obtener más información, consulte Arreglos altos.

Generación de código C/C++

Genere código C y C++ mediante MATLAB® Coder™.

Notas y limitaciones de uso:

• Los argumentos especificados utilizando argumentos nombre-valor deben ser constantes en tiempo de compilación.

• Los valores window de tamaño variable deben ser de doble precisión.

Entorno basado en subprocesos

Ejecute código en segundo plano con MATLAB® backgroundPool o acelere código con Parallel Computing Toolbox™ ThreadPool.

Notas y limitaciones de uso:

• La sintaxis sin argumentos de salida no es compatible.

• El vector de frecuencia debe estar espaciado uniformemente.

Para obtener más información, consulte Ejecutar funciones de MATLAB en un entorno basado en subprocesos.

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

R2023a: Visualizar salidas de funciones empleando la tarea de Live Editor Create Plot

Ahora puede usar la tarea de Live Editor Create Plot para visualizar la salida de spectrogram de forma interactiva. Puede seleccionar diferentes tipos de gráficas y establecer parámetros opcionales. La tarea también genera automáticamente código que se convierte en parte de su script en vivo.

Consulte también

Apps

• Signal Analyzer

Funciones

• goertzel | istft | periodogram | pspectrum | pwelch | stft | xspectrogram

Temas

• Spectrogram Computation with Signal Processing Toolbox
• Time-Frequency Gallery
• Practical Introduction to Time-Frequency Analysis
• Estimación de formantes con coeficientes LPC

Sitios web externos

• Análisis de Fourier (Recursos didácticos de MathWorks)

×

Abrir ejemplo

Tiene una versión modificada de este ejemplo. ¿Desea abrir este ejemplo con sus modificaciones?

No, sobrescribir la versión modificada Sí