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

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

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 точками ДПФ.

[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, путем lowpass фильтрации двух каналов и добавления их вместе. Задайте конечная импульсная характеристика 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 точками ДПФ. Постройте график оценки.

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 Гц. Одна из синусоид отстает от первого сигнала на Другая синусоида имеет задержку Оба сигнала встроены в белый Гауссов шум.

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 Гц. Таким образом, в третьей корреляционной функции нет peaks.

[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 Гц. Один из сигналов отстает друг от друга в фазе на  Встройте оба сигнала в белый Гауссов шум модуля отклонения.

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)