Main Content

spectrogram

Espectrograma utilizando la transformada de Fourier de tiempo corto

Descripción

ejemplo

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].

ejemplo

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

ejemplo

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

ejemplo

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

ejemplo

[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.

ejemplo

[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, [].

ejemplo

[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.

ejemplo

[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.

ejemplo

[___,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.

ejemplo

[___] = 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.

ejemplo

[___,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.

ejemplo

[___] = 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".

ejemplo

[___] = 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.

ejemplo

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

ejemplo

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

Ejemplos

contraer todo

Genere muestras Nx=1024 de una señal que consista en una suma de sinusoides. Las frecuencias normalizadas de las sinusoides son 2π/5 rad/muestra y 4π/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')

Figure contains an axes object. The axes object contains an object of type image.

Repita el cálculo.

  • Divida la señal en secciones de longitud nsc=Nx/4.5.

  • 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,2p) puntos, donde p=log2nsc.

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

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 M=49 muestras.

  • Especifique 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 NDFT=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)

Figure contains an axes object. The axes object contains an object of type patch.

Definición de STFT

Calcule la STFT de la señal de Nx muestras con la definición. Divida la señal en Nx-LM-L segmentos superpuestos. Aplique una ventana a cada segmento y evalúe la transformada discreta de Fourier en NDFT 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)

Figure contains an axes object. The axes object contains an object of type patch.

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

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.

subplot(2,1,1)
mesh(tp,fp,abs(sp).^2)
title("spectrogram")
view(2), axis tight
subplot(2,1,2)
mesh(t,f,abs(s).^2)
title("stft")
view(2), axis tight

Figure contains 2 axes objects. Axes object 1 with title spectrogram contains an object of type surface. Axes object 2 with title stft contains an object of type surface.

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 M=Nx/4.5, donde Nx 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(256,2log2M) 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 π en ambos casos. Para la salida de stft, divida los números de muestra por la tasa de muestreo efectiva, 2π.

subplot(2,1,1)
waterplot(sx,fx/pi,tx)
title("spectrogram")
subplot(2,1,2)
waterplot(st,ft/pi,tt/(2*pi))
title("stft")

Figure contains 2 axes objects. Axes object 1 with title spectrogram contains an object of type patch. Axes object 2 with title stft contains an object of type patch.

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

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 100/2+1=51.

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

Figure contains an axes object. The axes object contains an object of type image.

Realice un seguimiento de la frecuencia del chirp buscando la cresta de tiempo-frecuencia con más energía en los puntos de tiempo (2000-80)/(100-80)=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

Figure contains an axes object. The axes object contains 2 objects of type image, line.

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')

Figure contains an axes object. The axes object contains an object of type image.

Obtenga el mismo resultado calculando el espectrograma en 64 frecuencias equiespaciadas a lo largo del intervalo (-π,π]. 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')

Figure contains an axes object. The axes object contains an object of type image.

Genere una señal que conste de un oscilador controlado por tensión y tres átomos gaussianos. La señal se muestrea a fs=2 kHz durante 1 segundo.

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 Nx muestras en segmentos con longitud M=80 muestras, correspondientes a una resolución temporal de 80/2000=40 milisegundos.

  • Especifique 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 =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 g(n). El manchado y el factor de forma β de la ventana están relacionados por β=40×(1-).

  • pspectrum siempre usa NDFT=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 k=Nx-LM-L 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 ng(n) 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")

Figure contains 2 axes objects. Axes object 1 with title pspectrum contains an object of type patch. Axes object 2 with title spectrogram contains an object of type patch.

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 ng(n). 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)

Figure contains 2 axes objects. Axes object 1 with title pspectrum contains an object of type image. Axes object 2 with title spectrogram contains an object of type image.

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

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 β=18.

  • Especifique 120 muestras de solapamiento entre secciones contiguas.

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

spectrogram(y,kaiser(128,18),120,128,Fs,'reassigned','yaxis')

Figure contains an axes object. The axes object contains an object of type image.

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 β=18.

  • Especifique 120 muestras de solapamiento entre secciones contiguas.

  • Evalúe el espectro en frecuencias 128/2=65 y bins de tiempo (length(x)-120)/(128-120)=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),'.')

Figure contains an axes object. The axes object contains an object of type line.

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 M=73 muestras.

  • Especifique 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 NDFT=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");

Figure contains an axes object. The axes object contains an object of type image.

Calcule el espectrograma bilateral con la definición. Divida la señal en segmentos de M muestras con L muestras de solapamiento entre segmentos contiguos. Aplique una ventana a cada segmento y calcule la transformada discreta de Fourier en NDFT 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

subplot(2,1,1)
waterplot(Xtwo,ftwo,ttwo)
title("Two-Sided, Definition")
subplot(2,1,2)
waterplot(stwo,f,t)
title("Two-Sided, spectrogram Function")

Figure contains 2 axes objects. Axes object 1 with title Two-Sided, Definition contains an object of type patch. Axes object 2 with title Two-Sided, spectrogram Function contains an object of type patch.

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 NDFT, el intervalo de frecuencia está abierto en ambos extremos. En los valores pares de NDFT, 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

clf
subplot(2,1,1)
waterplot(Xcen,fcen,tcen)
title("Centered, Definition")
subplot(2,1,2)
waterplot(scen,f,t)
title("Centered, spectrogram Function")

Figure contains 2 axes objects. Axes object 1 with title Centered, Definition contains an object of type patch. Axes object 2 with title Centered, spectrogram Function contains an object of type patch.

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 NDFT, la STFT unilateral consiste en las primeras (NDFT+1)/2 filas de la STFT bilateral. En los valores pares de NDFT, la STFT unilateral consiste en las primeras NDFT/2+1 filas de la STFT bilateral.

  • En los valores impares de NDFT, el intervalo de frecuencia se cierra en la frecuencia cero y se abre en la frecuencia de Nyquist. En los valores pares de NDFT, 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

clf
subplot(2,1,1)
waterplot(Xone,fone,tone)
title("One-Sided, Definition")
subplot(2,1,2)
waterplot(sone,f,t)
title("One-Sided, spectrogram Function")

Figure contains 2 axes objects. Axes object 1 with title One-Sided, Definition contains an object of type patch. Axes object 2 with title One-Sided, spectrogram Function contains an object of type patch.

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

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])

Figure contains an axes object. The axes object contains an object of type patch.

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π.

[~,~,~,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

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

Figure contains an axes object. The axes object contains an object of type image.

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

Figure contains an axes object. The axes object contains 3 objects of type image, line.

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 β=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')

Figure contains an axes object. The axes object contains an object of type image.

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

view(-45,65)
colormap bone

Figure contains an axes object. The axes object contains an object of type surface.

Argumentos de entrada

contraer todo

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:

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.

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

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

  • Nw = window si window es un escalar.

  • Nw = length(window) si window es un 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]

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.

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.

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.

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.

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.

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

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

contraer todo

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 Nx, s tiene columnas k, donde

    • k = ⌊(Nxnoverlap)/(windownoverlap)⌋ si window es un escalar.

    • k = ⌊(Nxnoverlap)/(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".

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

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

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

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.

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

contraer todo

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 = MR 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=NxLML

 columnas, donde Nx 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 NDFT, 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 NDFT/2.

La columna m de la matriz de la STFT X(f)=[X1(f)X2(f)X3(f)Xk(f)] contiene la DFT de los datos con ventana centrados en el tiempo mR:

Xm(f)=n=x(n)g(nmR)ej2π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

Historial de versiones

Introducido antes de R2006a