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

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

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 выборок перекрытия между смежными сегментами и точками ДПФ 20:48.

[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. Окно сигналы с окном Hann с 1024 выборками. Задайте 512 выборок перекрытия между смежными сегментами и точками ДПФ 10:24. Постройте оценку.

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 Гц. Таким образом в третьей функции корреляции нет никакого 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 Гц. Один из сигналов изолирует другой в фазе π/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)