cpsd

Взаимная спектральная плотность мощности

Описание

пример

pxy = cpsd(x,y) оценивает взаимную спектральную плотность мощности (CPSD) двух сигналов дискретного времени, x и y, использование усредненного, модифицированного метода периодограммы валлийцев спектральной оценки.

  • Если x и y оба векторы, у них должна быть та же длина.

  • Если один из сигналов является матрицей, и другой вектор, то длина вектора должна равняться количеству строк в матрице. Функция расширяет вектор и возвращает матрицу оценок взаимной спектральной плотности мощности столбца столбцом.

  • Если x и y матрицы с одинаковым числом строк, но различными количествами столбцов, затем cpsd возвращает 3D массив, pxy, содержание взаимной спектральной плотности мощности оценивает для всех комбинаций входных столбцов. Каждый столбец pxy соответствует столбцу x, и каждая страница соответствует столбцу y: pxy(:,m,n) = cpsd(x(:,m),y(:,n)).

  • Если x и y матрицы равного размера, затем cpsd действует по столбцам: pxy(:,n) = cpsd(x(:,n),y(:,n)). Чтобы получить массив мультивхода/мультивыхода, добавьте 'mimo' к списку аргументов.

Для действительного x и y, cpsd возвращает односторонний CPSD. Для комплексного x или y, cpsd возвращает двухсторонний CPSD.

pxy = cpsd(x,y,window) использование window разделить x и y в сегменты и выполняют работу с окнами.

pxy = cpsd(x,y,window,noverlap) использование noverlap выборки перекрытия между смежными сегментами.

pxy = cpsd(x,y,window,noverlap,nfft) использование nfft выборка точек, чтобы вычислить дискретное преобразование Фурье.

пример

pxy = cpsd(___,'mimo') вычисляет массив мультивхода/мультивыхода оценок взаимной спектральной плотности мощности. Этот синтаксис может включать любую комбинацию входных параметров от предыдущих синтаксисов.

[pxy,w] = cpsd(___) возвращает вектор нормированных частот, w, в котором оценивается взаимная спектральная плотность мощности.

[pxy,f] = cpsd(___,fs) возвращает вектор частот, f, выраженный в терминах частоты дискретизации, fs, в котором оценивается взаимная спектральная плотность мощности. fs должен быть шестой числовой вход к cpsd. Чтобы ввести частоту дискретизации и все еще использовать значения по умолчанию предыдущих дополнительных аргументов, задайте эти аргументы как пустые, [].

пример

[pxy,w] = cpsd(x,y,window,noverlap,w) возвращает оценки взаимной спектральной плотности мощности на нормированных частотах, заданных в w.

пример

[pxy,f] = cpsd(x,y,window,noverlap,f,fs) возвращает оценки взаимной спектральной плотности мощности на частотах, заданных в f.

[___] = cpsd(x,y,___,freqrange) возвращает оценку взаимной спектральной плотности мощности по частотному диапазону, заданному freqrange. Допустимые опции для freqrange 'onesided', 'twosided', и 'centered'.

пример

cpsd(___) без выходных аргументов строит оценку взаимной спектральной плотности мощности в окне текущей фигуры.

Примеры

свернуть все

Сгенерируйте два цветных шумовых сигнала и постройте их взаимную спектральную плотность мощности. Задайте длину 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- демонстрационные экспоненциальные последовательности, xa=an и xb=bn, с n0. Задать 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;

Вычислите и постройте взаимную спектральную плотность мощности последовательностей на полном интервале нормированных частот, [-π,π]. Задайте прямоугольное окно длины N и никакое перекрытие между сегментами.

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

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

Перекрестный спектр мощности этих двух последовательностей имеет аналитическое выражение для большого N:

R(ω)=11-ae-jω11-bejω.

Преобразуйте это выражение во взаимную спектральную плотность мощности путем деления его на 2π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=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)

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

Входные параметры

свернуть все

Входные сигналы, заданные как векторы или матрицы.

Пример: cos(pi/4*(0:159))+randn(1,160) задает синусоиду, встроенную в белый Гауссов шум.

Типы данных: single | double
Поддержка комплексного числа: Да

Окно, заданное как целое число или как строка или вектор-столбец. Используйте window разделить сигнал на сегменты.

  • Если window целое число, затем cpsd делит x и y в сегменты длины window и окна каждый сегмент с Окном Хэмминга той длины.

  • Если window вектор, затем cpsd делит x и y в сегменты той же длины как вектор и окна каждый сегмент с помощью window.

Если длина x и y не может быть разделен точно в целое число сегментов с noverlap перекрывающиеся выборки, затем сигналы являются усеченными соответственно.

Если вы задаете window как пустой, затем cpsd использует Окно Хэмминга, таким образом что x и y разделены на восемь сегментов с noverlap перекрывающиеся выборки.

Для списка доступных окон смотрите Windows.

Пример: hann(N+1) и (1-cos(2*pi*(0:N)'/N))/2 оба задают окно Hann длины N  + 1.

Типы данных: single | double

Количество перекрытых выборок, заданных как положительное целое число.

  • Если window скаляр, затем noverlap должно быть меньшим, чем window.

  • Если window вектор, затем noverlap должно быть меньшим, чем длина window.

Если вы задаете noverlap как пустой, затем cpsd использует номер, который производит 50%-е перекрытие между сегментами. Если длина сегмента не задана, функция устанавливает noverlap к ⌊N/4.5 ⌋, где N является длиной сигналов ввода и вывода.

Типы данных: double | single

Количество точек ДПФ, заданных как положительное целое число. Если вы задаете nfft как пустой, затем cpsd устанавливает параметр на макс. (256,2p), где p = ⌈log2  N  для входных сигналов длины N.

Типы данных: single | double

Частота дискретизации, заданная как положительная скалярная величина. Частота дискретизации является количеством выборок в единицу времени. Если модуль времени является секундами, то частота дискретизации имеет модули Гц.

Нормированные частоты, заданные как строка или вектор-столбец по крайней мере с двумя элементами. Нормированные частоты находятся в рад/выборке.

Пример: w = [pi/4 pi/2]

Типы данных: double

Частоты, заданные как строка или вектор-столбец по крайней мере с двумя элементами. Частоты находятся в циклах в единицу времени. Единица времени задана частотой дискретизации, fs. Если fs имеет модули выборок/секунда, затем f имеет модули Гц.

Пример: fs = 1000; f = [100 200]

Типы данных: double

Частотный диапазон для оценки взаимной спектральной плотности мощности, заданной как 'onesided', 'twosided', или 'centered'. Значением по умолчанию является 'onesided' для сигналов с действительным знаком и 'twosided' для сигналов с комплексным знаком.

  • 'onesided' — Возвращает одностороннюю оценку взаимной спектральной плотности мощности двух входных сигналов с действительным знаком, x и y. Если nfft даже, pxy имеет nfft/2 + 1 строка и вычисляется на интервале [0, π] рад/выборка. Если nfft нечетно, pxy имеет (nfft + 1) строки/2 и интервал [0, π), рад/выборка. Если вы задаете fs, соответствующие интервалы [0, fs/2] циклы/единица времени для даже nfft и [0, fs/2) циклы/единица времени для нечетного nfft.

  • 'twosided' — Возвращает двухстороннюю оценку взаимной спектральной плотности мощности двух входных сигналов с комплексным знаком или с действительным знаком, x и y. В этом случае, pxy имеет nfft строки и вычисляются на интервале [0,2π), рад/выборка. Если вы задаете fs, интервал [0, fs) циклы/единица времени.

  • 'centered' — Возвращает двухстороннюю оценку в центре взаимной спектральной плотности мощности двух входных сигналов с комплексным знаком или с действительным знаком, x и y. В этом случае, pxy имеет nfft строки и вычисляются на интервале (–π, π] рад/выборка для даже nfft и (–π, π) рад/выборка для нечетного nfft. Если вы задаете fs, соответствующие интервалы (–fs/2, fs/2] циклы/единица времени для даже nfft и (–fs/2, fs/2) циклы/единица времени для нечетного nfft.

Выходные аргументы

свернуть все

Взаимная спектральная плотность мощности, возвращенная как вектор, матрица или 3D массив.

Нормированные частоты, возвращенные как вектор-столбец с действительным знаком.

  • Если pxy является односторонним, w охватывает интервал [0, π] когда nfft даже и [0, π), когда nfft нечетно.

  • Если pxy является двухсторонним, w охватывает интервал [0,2π).

  • Если pxy сосредоточен DC, w охватывает интервал (–π, π], когда nfft даже и (–π, π) когда nfft нечетно.

Типы данных: double | single

Частоты, возвращенные как вектор-столбец с действительным знаком.

Типы данных: double | single

Больше о

свернуть все

Взаимная спектральная плотность мощности

Взаимная спектральная плотность мощности является распределением степени на модульную частоту и задана как

Pxy(ω)=m=Rxy(m)ejωm.

Последовательность взаимной корреляции задана как

Rxy(m)=E{xn+myn}=E{xnynm},

где xn и yn являются совместно стационарными вероятностными процессами, – ∞ < n < ∞, <n<, и E {\·} оператор ожидаемого значения.

Алгоритмы

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

Ссылки

[1] Rabiner, Лоуренс Р. и Б. Голд. Теория и Приложение Цифровой обработки сигналов. Englewood Cliffs, NJ: Prentice Hall, 1975, стр 414–419.

[2] Валлийцы, Питер Д. “Использование Быстрого преобразования Фурье для Оценки Спектров мощности: Метод На основе Усреднения во времени По Коротким, Модифицированным Периодограммам”. IEEE® Transactions на Аудио и Электроакустике, издании AU-15, июнь 1967, стр 70–73.

[3] Оппенхейм, Алан V, Рональд В. Шафер и Джон Р. Бак. Обработка сигналов дискретного времени. 2-й Эд. Верхний Сэддл-Ривер, NJ: Prentice Hall, 1999.

Расширенные возможности

Генерация кода C/C++
Генерация кода C и C++ с помощью MATLAB® Coder™.

Представлено до R2006a

Для просмотра документации необходимо авторизоваться на сайте