Сгенерируйте два цветных шумовых сигнала и постройте их взаимную спектральную плотность мощности. Задайте длину 1 024 БПФ и треугольное окно с 500 точками без перекрытия.

r = randn(16384,1);

hx = fir1(30,0.2,rectwin(31));
x = filter(hx,1,r);

hy = ones(1,10)/sqrt(10);
y = filter(hy,1,r);

cpsd(x,y,triang(500),250,1024)

Сгенерируйте две двухканальных синусоиды, произведенные на уровне 1 кГц в течение 1 секунды. Каналы первого сигнала имеют частоты 200 Гц и 300 Гц. Каналы второго сигнала имеют частоты 300 Гц и 400 Гц. Оба сигнала встраиваются в белый Гауссов шум модульного отклонения.

fs = 1e3;
t = (0:1/fs:1-1/fs)';

q = 2*sin(2*pi*[200 300].*t);
q = q+randn(size(q));

r = 2*sin(2*pi*[300 400].*t);
r = r+randn(size(r));

Вычислите взаимную спектральную плотность мощности двух сигналов. Используйте окно Бартлетта с 256 выборками, чтобы разделить сигналы на сегменты и окно сегменты. Задайте 128 выборок перекрытия между смежными сегментами и точками ДПФ 20:48. Используйте встроенную функциональность cpsd построить результат.

cpsd(q,r,bartlett(256),128,2048,fs)

По умолчанию, cpsd работает столбец столбцом над матричными входными параметрами, одного размера. Каждый канал достигает максимума на частотах исходных синусоид.

Повторите вычисление, но теперь добавьте 'mimo' к списку аргументов.

cpsd(q,r,bartlett(256),128,2048,fs,'mimo')

Когда названо 'mimo' опция, cpsd возвращает 3D массив, содержащий оценки взаимной спектральной плотности мощности для всех комбинаций входных столбцов. Оценка второго канала q и первый канал r показывает расширенный пик на общей частоте 300 Гц.

Сгенерируйте два синусоидальных сигнала на 100 Гц, произведенные на уровне 1 кГц для 296 мс. Одна из синусоид изолирует другой на 2,5 мс, эквивалентные задержке фазы π/2. Оба сигнала встраиваются в белый Гауссов шум отклонения 1/4 ².

Fs = 1000;
t = 0:1/Fs:0.296;

x = cos(2*pi*t*100)+0.25*randn(size(t));
tau = 1/400;
y = cos(2*pi*100*(t-tau))+0.25*randn(size(t));

Вычислите и постройте величину взаимной спектральной плотности мощности. Используйте настройки по умолчанию для cpsd. Величина достигает максимума на частоте, где существует значительная когерентность между сигналами.

cpsd(x,y,[],[],[],Fs)

Постройте функцию квадрата когерентности и фазу перекрестного спектра. Ордината на частоте высокой когерентности соответствует задержке фазы между синусоидами.

[Cxy,F] = mscohere(x,y,[],[],[],Fs);

[Pxy,F] = cpsd(x,y,[],[],[],Fs);

subplot(2,1,1)
plot(F,Cxy)
title('Magnitude-Squared Coherence')

subplot(2,1,2)
plot(F,angle(Pxy))

hold on
plot(F,2*pi*100*tau*ones(size(F)),'--')
hold off

xlabel('Hz')
ylabel('\Theta(f)')
title('Cross Spectrum Phase')

Сгенерируйте два $$N$$- демонстрационные экспоненциальные последовательности, $$x_a=a^n$$ и $$x_b=b^n$$, с $$n\ge 0$$. Задать $$a=0.8$$, $$b=0.9$$, и маленькое $$N$$ видеть конечные размерные эффекты.

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

a = 0.8;
b = 0.9;

xa = a.^n;
xb = b.^n;

Вычислите и постройте взаимную спектральную плотность мощности последовательностей на полном интервале нормированных частот, $$[-\pi ,\pi ]$$. Задайте прямоугольное окно длины $$N$$ и никакое перекрытие между сегментами.

w = -pi:1/1000:pi;
wind = rectwin(N);
nove = 0;

[pxx,f] = cpsd(xa,xb,wind,nove,w);

Перекрестный спектр мощности этих двух последовательностей имеет аналитическое выражение для большого $$N$$:

Преобразуйте это выражение во взаимную спектральную плотность мощности путем деления его на $$2\pi N$$Сравнение результатов. Пульсация в cpsd результатом является последствие работы с окнами.

nfac = 2*pi*N;

X = 1./(1-a*exp(-1j*w));
Y = 1./(1-b*exp( 1j*w));
R = X.*Y/nfac;

semilogy(f/pi,abs(pxx))
hold on
semilogy(w/pi,abs(R))
hold off
legend('cpsd','Analytic')

Повторите вычисление с $$N=25$$. Кривые соглашаются на шесть фигур для $$N$$ всего 100.

N = 25;
n = 0:N-1;
xa = a.^n;
xb = b.^n;

wind = rectwin(N);

[pxx,f] = cpsd(xa,xb,wind,nove,w);
R = X.*Y/(2*pi*N);

semilogy(f/pi,abs(pxx))
hold on
semilogy(w/pi,abs(R))
hold off
legend('cpsd','Analytic')

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

Звуковые сигналы, сгенерированные, когда вы набираете номер или символ по цифровому телефону, являются суммами синусоид с частотами, взятыми из двух различных групп. Каждая пара тонов содержит одну частоту низкой группы (697 Гц, 770 Гц, 852 Гц или 941 Гц) и одну частоту высокой группы (1 209 Гц, 1 336 Гц или 1 477 Гц).

Сгенерируйте сигналы, соответствующие всем символам. Выборка каждый тон на уровне 4 кГц в течение половины секунды. Подготовьте таблицу ссылок.

fs = 4e3;
t = 0:1/fs:0.5-1/fs;

nms = ['1';'2';'3';'4';'5';'6';'7';'8';'9';'*';'0';'#'];

ver = [697 770 852 941];
hor = [1209 1336 1477];

v = length(ver);
h = length(hor);

for k = 1:v
    for l = 1:h
        idx = h*(k-1)+l;
        tone = sum(sin(2*pi*[ver(k);hor(l)].*t))';
        tones(:,idx) = tone;
    end
end

Постройте валлийскую периодограмму каждого сигнала и аннотируйте частоты компонента. Используйте Окно Хэмминга с 200 выборками, чтобы разделить сигналы на неперекрывающиеся сегменты и окно сегменты.

[pxx,f] = pwelch(tones,hamming(200),0,[],fs);

for k = 1:v
    for l = 1:h
        idx = h*(k-1)+l;
        ax = subplot(v,h,idx);
        plot(f,10*log10(pxx(:,idx)))
        ylim([-80 0])
        title(nms(idx))
        tx = [ver(k);hor(l)];
        ax.XTick = tx;
        ax.XTickLabel = int2str(tx);
    end
end

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

mys = sum(sin(2*pi*[ver(3);hor(2)].*t))'+5*randn(size(t'));

% To hear, type soundsc(mys,fs)

Вычислите взаимную спектральную плотность мощности поврежденного сигнала и опорных сигналов. Окно сигналы с помощью окна Кайзера с 512 выборками с масштабным фактором β = 5. Постройте величину каждого спектра.

[pxy,f] = cpsd(mys,tones,kaiser(512,5),100,[],fs);

for k = 1:v
    for l = 1:h
        idx = h*(k-1)+l;
        ax = subplot(v,h,idx);
        plot(f,10*log10(abs(pxy(:,idx))))
        ylim([-80 0])
        title(nms(idx))
        tx = [ver(k);hor(l)];
        ax.XTick = tx;
        ax.XTickLabel = int2str(tx);
    end
end

Цифра в поврежденном сигнале имеет спектр с самым высоким peaks и самым высоким значением RMS.

[~,loc] = max(rms(abs(pxy)));

digit = nms(loc)

digit = 
'8'