Процессы закона о степени

Важный класс сигналов с свободной от шкалы динамикой имеет автокорреляцию или степень спектральные плотности (PSD), которые следуют степени закону. Процесс закона о власти имеет PSD вида $$C|\omega {|}^{-\alpha }$$ для некоторой положительной константы, $$C$$, и некоторая экспонента $$\alpha$$. В некоторых случаях интересующий сигнал показывает PSD закона мощности. В других случаях интересующий сигнал повреждается шумом с PSD закона мощности. Эти шумы часто называют цветными. Умение оценить экспоненту от реализации этих процессов имеет важные последствия. Для одного, это позволяет вам сделать выводы о механизме, генерирующем данные, а также предоставить эмпирические доказательства для поддержки или отклонения теоретических предсказаний. В случае мешающего шума с PSD закона мощности, это полезно при разработке эффективных фильтров.

Коричневый шум, или брауновский процесс, является одним из таких цветных шумовых процессов с теоретическим показателем $$\alpha =2$$. Один из способов оценки экспоненты процесса закона степени - подгонка линии наименьших квадратов к логарифмическому графику PSD.

load brownnoise;
[Pxx,F] = pwelch(brownnoise,kaiser(1000,10),500,1000,1);
plot(log10(F(2:end)),log10(Pxx(2:end)));
grid on;
xlabel('log10(F)'); ylabel('log10(Pxx)');
title('Log-Log Plot of PSD Estimate')

Выполните регрессию значений PSD в журнале частот. Обратите внимание, что вы должны игнорировать нулевую частоту, чтобы избежать взятия нулевого логарифмического нуля.

Xpred = [ones(length(F(2:end)),1) log10(F(2:end))];
b = lscov(Xpred,log10(Pxx(2:end)));
y = b(1)+b(2)*log10(F(2:end));
hold on;
plot(log10(F(2:end)),y,'r--');
title(['Estimated Slope is ' num2str(b(2))]);

Кроме того, можно использовать как дискретный, так и непрерывный вейвлет методов анализа для оценки экспоненты. Отношения между экспонентом Держателя, H, возвращенные dwtleader и wtmm и $$\alpha$$ является ли этот сценарий $$\alpha =2H+1$$.

[dhbrown,hbrown,cpbrown] = dwtleader(brownnoise);
hexp = wtmm(brownnoise);
fprintf('Wavelet leader estimate is %1.2f\n',-2*cpbrown(1)-1);

Wavelet leader estimate is -1.91

fprintf('WTMM estimate is %1.2f\n',-2*hexp-1);

WTMM estimate is -2.00

В этом случае оценка, полученная путем подгонки линии наименьших квадратов к журналу оценки PSD, и оценки, полученной с помощью вейвлета методов, хорошо согласуются.

Мультифрактальный анализ

Существует ряд реальных сигналов, которые показывают нелинейное поведение закона мощности, которое зависит от моментов более высокого порядка и шкалы. Мультифрактальный анализ предоставляет способ описать эти сигналы. Многофракторный анализ состоит из определения, существует ли какой-либо тип масштабирования закона мощности для различных статистических моментов в разных шкалах. Если это поведение масштабирования характеризуется одной степенью масштабирования или эквивалентно является линейной функцией моментов, процесс является монофрактальным. Если поведение масштабирования по шкале является нелинейной функцией моментов, процесс является многофракторным. Коричневый шум из предыдущего раздела является примером монофрактального процесса, и это показано в более позднем разделе.

Чтобы проиллюстрировать, как фрактальный анализ может выявить структуру сигнала, не очевидную с помощью более классических методов обработки сигнала, загружайте RWdata.mat который содержит два временных рядов (Ts1 и Ts2) с 8000 выборок каждый. Постройте график данных.

load RWdata;
figure;
plot([Ts1 Ts2]); grid on;
legend('Ts1','Ts2','Location','NorthEast');
xlabel('Time'); ylabel('Amplitude');

Сигналы имеют очень похожую статистику второго порядка. Если вы посмотрите на средства, значения RMS и отклонения Ts1 и Ts2, значения почти идентичны. Оценки PSD также очень похожи.

pwelch([Ts1 Ts2],kaiser(1e3,10))

Автокорреляционные последовательности распадаются очень медленно для обоих временных рядов и не информативны для дифференцирования временных рядов.

[xc1,lags] = xcorr(Ts1,300,'coef');
xc2 = xcorr(Ts2,300,'coef');
subplot(2,1,1)
hs1 = stem(lags(301:end),xc1(301:end));
hs1.Marker = 'none';
title('Autocorrelation Sequence of Ts1');
subplot(2,1,2)
hs2 = stem(lags(301:end),xc2(301:end));
hs2.Marker = 'none';
title('Autocorrelation Sequence of Ts2');
xlabel('Lag')

Даже при задержке 300 автокорреляции составляют 0,94 и 0,96 соответственно.

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

[dh1,h1,cp1,tauq1] = dwtleader(Ts1);
[dh2,h2,cp2,tauq2] = dwtleader(Ts2);
figure;
hp = plot(h1,dh1,'b-o',h2,dh2,'b-^');
hp(1).MarkerFaceColor = 'b';
hp(2).MarkerFaceColor = 'r';
grid on;
xlabel('h'); ylabel('D(h)');
legend('Ts1','Ts2','Location','NorthEast');
title('Multifractal Spectrum');

Мультифрактальный спектр эффективно показывает распределение показателей масштабирования для сигнала. Эквивалентно, мультифрактальный спектр обеспечивает меру того, насколько локальная регулярность сигнала изменяется во времени. Сигнал, который является монофрактальным, проявляет по существу одинаковую регулярность повсюду во времени и, следовательно, имеет мультифрактальный спектр с узкой поддержкой. И наоборот, мультифрактальный сигнал показывает изменения регулярности сигнала с течением времени и имеет мультифрактальный спектр с более широкой поддержкой. Из показанных здесь мультифрактальных спектров, Ts2, по-видимому, является монофрактальным сигналом, характеризующимся кластером масштабирующих показателей около 0,78. С другой стороны Ts1, демонстрирует широкий диапазон масштабирующих показателей, указывающих, что он является многофракторным. Обратите внимание на общую область значений показателей масштабирования (Holder) для Ts2 всего 0,14, в то время как это в 4,6 раза больше, чем для Ts1. Ts2 является на самом деле примером монофрактального дробного процесса Броуновского движения (fBm) с степенью Холдера 0,8 и Ts1 является многофракционной случайной ходьбой.

Можно также использовать выходные параметры масштабирования экспоненты от dwtleader наряду со 2-м кумулянтом, чтобы помочь классифицировать процесс как монофрактальный по сравнению с мультифрактальным. Напомним, что монофрактальный процесс имеет линейный закон масштабирования как функцию от статистических моментов, в то время как мультифрактальный процесс имеет нелинейный закон масштабирования. dwtleader использует область значений моментов от -5 до 5 при оценке этих законов масштабирования. График показателей масштабирования для процесса fBm и многофракционной случайной ходьбы (MRW) показывает различие.

plot(-5:5,tauq1,'b-o',-5:5,tauq2,'r-^');
grid on;
xlabel('Q-th Moment'); ylabel('Scaling Exponents');
title('Scaling Exponents');
legend('MRW','fBm','Location','SouthEast');

Экспоненты масштабирования для процесса fBm являются линейной функцией моментов, в то время как экспоненты для процесса MRW показывают отклонение от линейности. Ту же информацию суммируют 1-я, 2-я и 3-я кумулянты. Первый кумулянт является оценкой наклона, другими словами, он захватывает линейное поведение. Второй кумулянт захватывает первый уклон от линейности. Можно думать о втором кумулянте как о коэффициентах для члена второго порядка (квадратичного), в то время как третий кумулянт характеризует более сложное отклонение масштабирующих экспонентов от линейности. Если вы исследуете 2-й и 3-й кумулянты для процесса MRW, они в 6 и 42 раза больше, чем соответствующие кумулянты для данных fBm. В последнем случае 2-й и 3-й кумулянты почти равны нулю, как и ожидалось.

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

hp = plot(h1,dh1,'b-o',h2,dh2,'b-^',hbrown,dhbrown,'r-v');
hp(1).MarkerFaceColor = 'b';
hp(2).MarkerFaceColor = 'r';
hp(3).MarkerFaceColor = 'k';
grid on;
xlabel('h'); ylabel('D(h)');
legend('Ts1','Ts2','brown noise','Location','SouthEast');
title('Multifractal Spectrum');

Где процесс идет дальше? Стойкое и антиперсистентное поведение

Оба дробных брауновских процесса (Ts2) и коричневый шумовой ряд являются монофрактальными. Однако простой график двух временных рядов показывает, что они появляются совсем по-другому.

subplot(2,1,1)
plot(brownnoise); title('Brown Noise');
grid on; axis tight;
subplot(2,1,2)
plot(Ts2); title('fBm H=0.8'); grid on;
axis tight;

Данные fBm намного плавнее, чем коричневый шум. Шум Брауна, также известный как случайная прогулка, имеет теоретический показатель Холдера 0,5. Это значение формирует контур между процессами с показателями Держателя, H, от 0 < H < 0,5 и процессами с показателями Держателя в интервале 0,5 < H < 1. Первые называются антиперсистентными и демонстрируют короткую память. Последние называются стойкими и демонстрируют долгую память. В антиперсистентных временных рядах следует увеличение значения в момент t с высокой вероятностью с уменьшением значения в момент t + 1. Точно так же уменьшение значения в момент t обычно сопровождается увеличением значения в момент t + 1. Другими словами, временные ряды имеет тенденцию всегда возвращаться к своему среднему значению. В постоянных временных рядах увеличения значения, как правило, сопровождаются последующими увеличениями, в то время как уменьшение значения, как правило, сопровождается последующим уменьшением.

Чтобы увидеть некоторые реальные примеры антиперсистентных временных рядов, загрузите и проанализируйте ежедневные журналы возвратов для фондовых индексов Тайваня и Сеула. Дневные возвраты по обоим индексам охватывают приблизительный период с июля 1997 года по апрель 2016 года.

load StockCompositeData;
subplot(2,1,1)
plot(SeoulComposite); title('Seoul Composite Index - 07/1997-04/2016');
ylabel('Log Returns'); grid on;
subplot(2,1,2);
plot(TaiwanWeighted); title('Taiwan Weighted Index - 07/1997-04/2016');
ylabel('Log Returns');
grid on;

Получите и постройте график многофракционных спектров этих двух временных рядов.

[dhseoul,hseoul,cpseoul] = dwtleader(SeoulComposite);
[dhtaiwan,htaiwan,cptaiwan] = dwtleader(TaiwanWeighted);
figure;
plot(hseoul,dhseoul,'b-o','MarkerFaceColor','b');
hold on;
plot(htaiwan,dhtaiwan,'r-^','MarkerFaceColor','r');
xlabel('h'); ylabel('D(h)'); grid on;
title('Multifractal Spectrum');

Из мультифрактального спектра ясно, что оба временных рядов являются антиперсистентными. Для сравнения постройте график многофракционных спектров двух финансовых временных рядов вместе с данными о коричневом шуме и fBm, показанными ранее.

plot(hbrown,dhbrown,'k-v','MarkerFaceColor','k');
plot(h2,dh2,'b-*','MarkerFaceColor','b');
legend('Seoul Composite','Taiwan Weighted Index','Brown Noise','FBM',...
    'Location','SouthEast');
hold off;

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

Измерение фрактальной динамики вариабельности сердечного ритма

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

В качестве примера индуцированного изменения фрактальной динамики изменчивости сердечного ритма рассмотрим пациент, вводимый E1 простагландина из-за тяжелого гипертонического эпизода. Данные являются частью RHRV, основанного на R программного пакета для анализа изменчивости сердечного ритма. Авторы любезно предоставили разрешение на его использование в этом примере.

Загрузите и постройте график данных. Вертикальная красная линия отмечает начало эффекта простагландина на E1 сердечного ритма и вариабельности сердечного ритма.

load hrvDrug;
plot(hrvDrug); grid on;
hold on;
plot([4642 4642],[min(hrvDrug) max(hrvDrug)],'r','linewidth',2);
hold off;
ylabel('Heart Rate'); xlabel('Sample');

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

predrug = hrvDrug(1:4642);
postdrug = hrvDrug(4643:end);
[dhpre,hpre] = dwtleader(predrug);
[dhpost,hpost] = dwtleader(postdrug);
figure;
hl = plot(hpre,dhpre,'b-d',hpost,dhpost,'r-^');
hl(1).MarkerFaceColor = 'b';
hl(2).MarkerFaceColor = 'r';
xlabel('h'); ylabel('D(h)');
grid on;
legend('Predrug','Postdrug');
title('Multifractal Spectrum'); xlabel('h'); ylabel('D(h)');

Индукция препарата привела к 50% снижению ширины фрактального спектра. Это указывает на значительное снижение нелинейной динамики сердца, измеренное изменчивостью сердечного ритма. В этом случае уменьшение фрактальной размерности было частью медицинского вмешательства. В другом контексте исследования групп здоровых индивидуумов и пациентов с застойным сердечным отказом показали, что различия в многофракционных спектрах могут дифференцировать эти группы. В частности, значительное снижение ширины мультифрактального спектра является маркером сердечной дисфункции.

Ссылки

L. Rodriguez-Linares, L., A.J. Mendez, M.J. Lado, D.N. Olivieri, X.A. Vila, and I. Gomez-Conde, «Открытый исходный код для спектрального анализа вариабельности сердечного ритма», Компьютерные методы и программы в биологии

Wendt, H. and Abry, P. «Multifractality test using bootstrapped wavelet leaders», IEEE Trans. Signal Processing, vol. 55, no. 10, pp. 4811-4820, 2007.

Wendt, H., Abry, P. and Jaffard, S. «Bootstrap for empirical multifractal analysis», IEEE Signal Processing Magazine, 24, 4, 38-48, 2007.

Jaffard, S., Lashermes, B., and Abry, P. «Wavelet leaders in multifractal analysis». В Т. Цянь, М. И. Ваи и X. Юэшэн, редактора. Анализ и применение вейвлет, стр. 219-264, Биркхаузер, 2006.