Генерация $${}_{\mathrm{Nx}}\mathrm{=}\mathrm{}\mathrm{}\mathrm{}$$1024 выборок сигнала, который состоит из суммы синусоид. Нормализованные частоты синусоид $$\mathrm{составляют}\mathrm{}$$2δ/5 рад/образец $$\mathrm{и}\mathrm{}$$4λ/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 точек DFT и окно Хэмминга по умолчанию.

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 точек DFT.

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 точек DFT.

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{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}$$100/2+1⌋=51 частотах.

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

Отслеживайте частоту чирпа, обнаруживая частотно-временной гребень с наибольшей энергией по $$\lfloor \mathrm{(}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}2000-80\mathrm{}\mathrm{)/(}\mathrm{}\mathrm{}\mathrm{}100-80\mathrm{)}\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 точки DFT. Постройте график спектрограммы.

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

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

Получим тот же результат вычислением спектрограммы на 64 равнопроходных частотах на интервале $(-\lambda ,\lambda$]. '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, которые меньше $$-\mathrm{}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 с параметром формы $$\beta \mathrm{=}\mathrm{}$$17, 10 выборок перекрытия между соседними секциями и длиной БПФ 256.

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

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

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

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

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

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

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

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 \mathrm{=}$$5. Укажите 220 образцов перекрытия сечения и 512 точек DFT. Постройте график частоты по оси Y. Используйте карту цветов и представление по умолчанию.

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

Измените вид для отображения спектрограммы в виде графика водопада. Задайте для карты цветов значение bone.

view(-45,65)
colormap bone