Вычислите и постройте график оценки когерентности между двумя цветными шумовыми последовательностями.

Генерировать сигнал, состоящий из белого гауссова шума.

r = randn(16384,1);

Для создания первой последовательности полосовой фильтр сигнала. Спроектируйте фильтр 16-го порядка, который пропускает нормализованные частоты между 0,2δ и 0,4δ рад/выборка. Задайте затухание полосы останова 60 дБ. Фильтрация исходного сигнала.

dx = designfilt('bandpassiir','FilterOrder',16, ...
    'StopbandFrequency1',0.2,'StopbandFrequency2',0.4, ...
    'StopbandAttenuation',60);
x = filter(dx,r);

Чтобы создать вторую последовательность, сконструируйте фильтр 16-го порядка, который останавливает нормализованные частоты между 0.6δ и 0.8δ рад/выборка. Укажите пульсацию полосы пропускания 0,1 дБ. Фильтрация исходного сигнала.

dy = designfilt('bandstopiir','FilterOrder',16, ...
    'PassbandFrequency1',0.6,'PassbandFrequency2',0.8, ...
    'PassbandRipple',0.1);
y = filter(dy,r);

Оценка квадратичной по величине когерентности x и y. Используйте 512-образное окно Хэмминга. Укажите 500 выборок перекрытия между соседними сегментами и 2048 точками DFT.

[cxy,fc] = mscohere(x,y,hamming(512),500,2048);

Постройте график функции когерентности и наложите частотные характеристики фильтров.

[qx,f] = freqz(dx);
qy = freqz(dy);

plot(fc/pi,cxy)
hold on
plot(f/pi,abs(qx),f/pi,abs(qy))
hold off

Генерировать случайный двухканальный сигнал, x. Генерировать другой сигнал, y, путем фильтрации нижних частот двух каналов и их сложения. Укажите фильтр FIR 30-го порядка с частотой отсечения 0,3δ, разработанный с использованием прямоугольного окна.

h = fir1(30,0.3,rectwin(31));
x = randn(16384,2);
y = sum(filter(h,1,x),2);

Вычислить оценку множественной когерентности x и y. Оконные сигналы с 1024-образным окном Ханна. Укажите 512 образцов перекрытия между соседними сегментами и 1024 точками DFT. Постройте график оценки.

noverlap = 512;
nfft = 1024;

mscohere(x,y,hann(nfft),noverlap,nfft,'mimo')

Сравните оценку когерентности с частотной характеристикой фильтра. Падения когерентности соответствуют нулям частотной характеристики.

[H,f] = freqz(h);

hold on
yyaxis right
plot(f/pi,20*log10(abs(H)))
hold off

Вычислите и постройте график обычной квадратичной оценки когерентности x и y. Оценка не достигает 1 ни для одного из каналов.

figure
mscohere(x,y,hann(nfft),noverlap,nfft)

Генерируйте два многоканальных сигнала, каждый из которых дискретизируется при частоте 1 кГц в течение 2 секунд. Первый сигнал, входной, состоит из трёх синусоид с частотами 120 Гц, 360 Гц и 480 Гц. Второй сигнал, выходной, состоит из двух синусоид с частотами 120 Гц и 360 Гц. Одна из синусоид задерживает первый сигнал на δ/2. Другая синусоида имеет отставание δ/4. Оба сигнала встроены в белый гауссов шум.

fs = 1000;
f = 120;
t = (0:1/fs:2-1/fs)';

inpt = sin(2*pi*f*[1 3 4].*t);
inpt = inpt+randn(size(inpt));
oupt = sin(2*pi*f*[1 3].*t-[pi/2 pi/4]);
oupt = oupt+randn(size(oupt));

Оцените степень корреляции между всеми входными сигналами и каждым из выходных каналов. Для отображения данных используйте окно Хэмминга длиной 100. mscohere возвращает одну функцию когерентности для каждого выходного канала. Функции когерентности достигают максимумов на частотах, совместно используемых входом и выходом.

[Cxy,f] = mscohere(inpt,oupt,hamming(100),[],[],fs,'mimo');

for k = 1:size(oupt,2)
    subplot(size(oupt,2),1,k)
    plot(f,Cxy(:,k))
    title(['Output ' int2str(k) ', All Inputs'])
end

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

[Cxy,f] = mscohere(oupt,inpt,hamming(100),[],[],fs,'mimo');

for k = 1:size(inpt,2)
    subplot(size(inpt,2),1,k)
    plot(f,Cxy(:,k))
    title(['Input ' int2str(k) ', All Outputs'])
end

Повторите вычисления с помощью функции печати mscohere.

clf
mscohere(oupt,inpt,hamming(100),[],[],fs,'mimo')

Вычисляют обычную функцию когерентности второго сигнала и первых двух каналов первого сигнала. Непиковые значения отличаются от функции множественной когерентности.

[Cxy,f] = mscohere(oupt,inpt(:,[1 2]),hamming(100),[],[],fs);
plot(f,Cxy)

Найдите разности фаз, вычисляя угол перекрестного спектра в точках максимальной когерентности.

Pxy = cpsd(oupt,inpt(:,[1 2]),hamming(100),[],[],fs);
[~,mxx] = max(Cxy);
for k = 1:2
    fprintf('Phase lag %d = %5.2f*pi\n',k,angle(Pxy(mxx(k),k))/pi)
end

Phase lag 1 = -0.51*pi
Phase lag 2 = -0.22*pi

Генерируют два синусоидальных сигнала, дискретизированных в течение 1 секунды каждый на частоте 1 кГц. Каждая синусоида имеет частоту 250 Гц. Один из сигналов отстает от другого по фазе на δ/3 радиан. Встроить оба сигнала в белый гауссов шум единичной дисперсии.

fs = 1000;
f = 250;
t = 0:1/fs:1-1/fs;
um = sin(2*pi*f*t)+rand(size(t));
un = sin(2*pi*f*t-pi/3)+rand(size(t));

Использовать mscohere для вычисления и построения графика квадратичной когерентности сигналов.

mscohere(um,un,[],[],[],fs)

Измените название графика, метку оси X и границы оси Y.

title('Magnitude-Squared Coherence')
xlabel('f (Hz)')
ylim([0 1.1])

Использовать gca для получения маркера перемещения к текущим осям. Измените расположение делений. Удалите метку оси Y.

ax = gca;
ax.XTick = 0:250:500;
ax.YTick = 0:0.25:1;
ax.YLabel.String = [];

Позвоните в Children свойство маркера перемещения для изменения цвета и ширины линии печати.

ln = ax.Children;
ln.Color = [0.8 0 0];
ln.LineWidth = 1.5;

В качестве альтернативы можно использовать set и get для изменения свойств линии.

set(get(gca,'Children'),'Color',[0 0.4 0],'LineStyle','--','LineWidth',1)