Сгенерируйте два цветных шумовых сигнала и постройте их перекрестную степень спектральная плотность. Задайте длину 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=\mathrm{}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)

Вычислите перекрестную степень спектральная плотность поврежденного сигнала и ссылочных сигналов. Окно сигналы с помощью окна Kaiser с 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'