Сгенерируйте два цветных сигнала шума и постройте график их взаимной спектральной плотности мощности. Задайте БПФ длины 1024 и треугольное окно с 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-выборочное окно Bartlett, чтобы разделить сигналы на сегменты и отобразить сегменты. Задайте 128 выборок перекрытия между смежными сегментами и 2048 точками ДПФ. Используйте встроенную функциональность cpsd для построения графика результатов.

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

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

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

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

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

Сгенерировать два синусоидальных сигнала 100 Гц, дискретизированных с частотой дискретизации 1 кГц в течение 296 мс. Оба сигнала встроены в белый Гауссов шум отклонения 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 Гц) и одну частоту высокой группы (1209 Гц, 1336 Гц или 1477 Гц).

Сгенерируйте сигналы, соответствующие всем символам. Дискретизируйте каждый тон с частотой 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

Постройте график периодограммы Welch каждого сигнала и аннотируйте частоты компонентов. Используйте 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'