Сгенерировать $$N_x=\mathrm{}\mathrm{}\mathrm{}1024$$ выборки сигнала, который состоит из суммы синусоид. Нормированные частоты синусоид $$2\pi /5$$ рад/выборка и $$4\pi /5$$ рад/выборка. Более высокая синусоида частоты имеет 10 раз амплитуду другой синусоиды.

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

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

Вычислите кратковременное преобразование Фурье с помощью функциональных значений по умолчанию. Постройте спектрограмму.

s = spectrogram(x);

spectrogram(x,'yaxis')

Повторите вычисление.

Проверьте, что два подхода дают идентичные результаты.

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

Разделите сигнал на 8 разделов равной длины с 50%-м перекрытием между разделами. Задайте ту же длину БПФ как на предыдущем шаге. Вычислите кратковременное преобразование Фурье и проверьте, что оно дает тот же результат как предыдущие две процедуры.

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

Сгенерируйте квадратичный щебет, x, выбранный на уровне 1 кГц в течение 2 секунд. Частота щебета составляет 100 Гц первоначально и пересекает 200 Гц в t = 1 с.

t = 0:0.001:2;
x = chirp(t,100,1,200,'quadratic');

Вычислите и отобразите спектрограмму x.

spectrogram(x,128,120,128,1e3)

Замените Окно Хэмминга на окно Блэкмена. Уменьшите перекрытие к 60 выборкам. Постройте ось времени так, чтобы ее значения увеличились сверху донизу.

spectrogram(x,blackman(128),60,128,1e3)
ax = gca;
ax.YDir = 'reverse';

Вычислите и отобразите PSD каждого сегмента квадратичного щебета, который запускается на уровне 100 Гц и пересекает 200 Гц в t = 1 секунда. Задайте частоту дискретизации 1 кГц, продолжительность сегмента 128 выборок и перекрытие 120 выборок. Используйте 128 точек ДПФ и Окно Хэмминга по умолчанию.

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

spectrogram(x,128,120,128,fs,'yaxis')
title('Quadratic Chirp')

Вычислите и отобразите PSD каждого сегмента линейного щебета, выбранного на уровне 1 кГц, который запускается в DC и пересекает 150 Гц в t = 1 секунда. Задайте продолжительность сегмента 256 выборок и перекрытие 250 выборок. Используйте Окно Хэмминга по умолчанию и 256 точек ДПФ.

x = chirp(t,0,1,150);

spectrogram(x,256,250,256,fs,'yaxis')
title('Linear Chirp')

Вычислите и отобразите PSD каждого сегмента логарифмического щебета, выбранного на уровне 1 кГц, который запускается на уровне 20 Гц и пересекает 60 Гц в t = 1 секунда. Задайте продолжительность сегмента 256 выборок и перекрытие 250 выборок. Используйте Окно Хэмминга по умолчанию и 256 точек ДПФ.

x = chirp(t,20,1,60,'logarithmic');

spectrogram(x,256,250,[],fs,'yaxis')
title('Logarithmic Chirp')

Используйте логарифмический масштаб для оси частоты. Спектрограмма становится строкой.

ax = gca;
ax.YScale = 'log';

Используйте функцию spectrogram, чтобы измерить и отследить мгновенную частоту сигнала.

Сгенерируйте квадратичный щебет, выбранный на уровне 1 кГц в течение двух секунд. Задайте щебет так, чтобы его частота была первоначально 100 Гц и увеличилась до 200 Гц после одной секунды.

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

Оцените спектр щебета с помощью кратковременного преобразования Фурье, реализованного в функции spectrogram. Разделите сигнал на разделы длины 100, оконный с Окном Хэмминга. Задайте 80 выборок перекрытия между смежными разделами и оцените спектр в $$\lfloor \mathrm{}\mathrm{}100/2+1\rfloor =\mathrm{}51$$ частоты.

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

Отследите частоту щебета путем нахождения гребня частоты времени с самой высокой энергией через $$\lfloor (\mathrm{}\mathrm{}\mathrm{}2000-\mathrm{}80)/(\mathrm{}\mathrm{}100-\mathrm{}80)\rfloor =\mathrm{}96$$ моменты времени. Наложите мгновенную частоту на графике спектрограммы.

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

Сгенерируйте 512 выборок щебета с синусоидально переменным содержимым частоты.

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

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

Вычислите двухстороннее кратковременное преобразование Фурье в центре щебета. Разделите сигнал на сегменты с 32 выборками с перекрытием с 16 выборками. Задайте 64 точки ДПФ. Постройте спектрограмму.

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

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

Получите тот же результат путем вычисления спектрограммы на 64 равномерно расположенных частотах на интервале $(-\pi ,\pi ]$. Опция 'centered' не необходима.

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

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

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

Сгенерируйте сигнал щебета, выбранный в течение 2 секунд на уровне 1 кГц. Задайте щебет так, чтобы его частота была первоначально 100 Гц и увеличилась до 200 Гц после 1 секунды.

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

Оцените переприсвоенную спектрограмму сигнала.

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

Сгенерируйте сигнал щебета, выбранный в течение 2 секунд на уровне 1 кГц. Задайте щебет так, чтобы его частота была первоначально 100 Гц и увеличилась до 200 Гц после 1 секунды.

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

Оцените зависящую от времени степень спектральную плотность (PSD) сигнала.

Выведите частоту и время центра тяжести каждой оценки PSD. Обнулите те элементы меньшего PSD, чем $$-\mathrm{}30$$ дБ.

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

Постройте ненулевые элементы как функции частот центра тяжести и времена.

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

Сгенерируйте сигнал, выбранный на уровне 1 024 Гц в течение 2 секунд.

nSamp = 2048;
Fs = 1024;
t = (0:nSamp-1)'/Fs;

В течение первой секунды сигнал состоит из синусоиды на 400 Гц и вогнутого квадратичного щебета. Задайте щебет так, чтобы это было симметрично о средней точке интервала, начав и закончившись на частоте 250 Гц и достигнув минимума 150 Гц.

t1 = t(1:nSamp/2);

x11 = sin(2*pi*400*t1);
x12 = chirp(t1-t1(nSamp/4),150,nSamp/Fs,1750,'quadratic');
x1 = x11+x12;

Остальная часть сигнала состоит из двух линейных щебетов уменьшающейся частоты. Один щебет имеет начальную частоту 250 Гц, которая уменьшается до 100 Гц. Другой щебет имеет начальную частоту 400 Гц, которая уменьшается до 250 Гц.

t2 = t(nSamp/2+1:nSamp);

x21 = chirp(t2,400,nSamp/Fs,100);
x22 = chirp(t2,550,nSamp/Fs,250);
x2 = x21+x22;

Добавьте белый Гауссов шум в сигнал. Задайте отношение сигнал-шум 20 дБ. Сбросьте генератор случайных чисел для восстанавливаемых результатов.

SNR = 20;
rng('default')

sig = [x1;x2];
sig = sig + randn(size(sig))*std(sig)/db2mag(SNR);

Вычислите и постройте спектрограмму сигнала. Задайте окно Kaiser длины 63 с параметром формы $$\beta =\mathrm{}17$$, 10 меньше выборок перекрытия между смежными разделами и длина БПФ 256.

nwin = 63;
wind = kaiser(nwin,17);
nlap = nwin-10;
nfft = 256;

spectrogram(sig,wind,nlap,nfft,Fs,'yaxis')

Порог спектрограмма так, чтобы были обнулены любые элементы со значениями, меньшими, чем ОСШ.

spectrogram(sig,wind,nlap,nfft,Fs,'MinThreshold',-SNR,'yaxis')

Повторно присвойте каждую оценку PSD местоположению ее центра энергии.

spectrogram(sig,wind,nlap,nfft,Fs,'reassign','yaxis')

Порог переприсвоенная спектрограмма так, чтобы были обнулены любые элементы со значениями, меньшими, чем ОСШ.

spectrogram(sig,wind,nlap,nfft,Fs,'reassign','MinThreshold',-SNR,'yaxis')

Загрузите звуковой сигнал, который содержит два уменьшающихся щебета, и широкополосное обрызгивают звук. Вычислите кратковременное преобразование Фурье. Разделите форму волны на сегменты с 400 выборками с перекрытием с 300 выборками. Постройте спектрограмму.

load splat

% To hear, type soundsc(y,Fs)

sg = 400;
ov = 300;

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

Используйте функцию spectrogram, чтобы вывести степень спектральную плотность (PSD) сигнала.

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

Отследите два щебета с помощью функции medfreq. Чтобы найти более сильный, низкочастотный щебет, ограничьте поиск частотами выше 100 Гц и временами перед запуском широкополосного звука.

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

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

Чтобы найти слабый высокочастотный щебет, ограничьте поиск частотами выше 2 500 Гц и временами между 0,3 секундами и 0,65 секундами.

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

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

Наложите результат на спектрограмме. Разделите значения частоты на 1 000, чтобы выразить их в kHz.

hold on
plot(t(t1),m1/1000,'linewidth',4)
plot(t(t2),m2/1000,'linewidth',4)
hold off

Сгенерируйте две секунды сигнала, выбранного на уровне 10 кГц. Задайте мгновенную частоту сигнала как треугольная функция времени.

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

Вычислите и постройте спектрограмму сигнала. Используйте окно Kaiser длины 256 и сформируйте параметр $$\beta =5$$. Задайте 220 выборок перекрытия от раздела к разделу и 512 точек ДПФ. Постройте частоту на оси Y. Используйте палитру по умолчанию и представление.

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

Измените представление, чтобы отобразить спектрограмму как график водопада. Установите палитру на bone.

colormap bone
view(-45,65)