Сгенерировать $$N_x=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 100/2+1\rfloor =51$$ частоты.

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

Отследите частоту щебета путем нахождения частотно-временного гребня с самой высокой энергией через $$\lfloor (2000-80)/(100-80)\rfloor =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, чем $$-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);

Вычислите и постройте спектрограмму сигнала. Задайте окно Кайзера длины 63 параметром формы $$\beta =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);

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

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

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

view(-45,65)
colormap bone