Произвести $$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 кГц, который начинается с постоянного тока и пересекает 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 ]$. The '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),'.')

Сгенерируйте сигнал, дискретизированный на частоте 1024 Гц в течение 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 с параметром shape $$\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));

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

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

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

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

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