Создайте два цветных шумовых сигнала и постройте график их спектральной плотности. Задайте длину - 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-образное окно Бартлетта, чтобы разделить сигналы на сегменты и открыть сегменты. Укажите 128 выборок перекрытия между соседними сегментами и 2048 точками DFT. Использовать встроенные функциональные возможности 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 мс. Одна из синусоид отстает от другой на 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-образцов, $${}_{\mathrm{xa}}{=}^{}$$an $${и}_{}{\mathrm{xb}}^{}$$= bn, $$\mathrm{с}$$n≥0. $$\mathrm{Задайте}\mathrm{a}$$= $$\mathrm{}0,8\mathrm{,}$$b = 0,9 и $$$$небольшое значение N, чтобы увидеть конечные эффекты.

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

a = 0.8;
b = 0.9;

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

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

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

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

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

Преобразуйте это выражение в спектральную плотность поперечной мощности, деля его на $$\mathrm{}\mathrm{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\mathrm{=}\mathrm{}$$25. Кривые согласуются с шестью цифрами $$\mathrm{для}$$ 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

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

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

digit = nms(loc)

digit = 
'8'