Анализируйте изменчивость сигнала по шкале

Загрузите и постройте график clock_571 набор данных. Этот пример по существу является воссозданием анализа, представленного в Percival & Walden [3], стр. 13-16.

load clock_571;
figure;
plot(clock_571)
xlabel('Days')
grid on;
title('Daily Average Fractional Frequency Deviates -- Cesium Clock');

Данные представляют собой ежедневные отклонения дробной частоты среднего значения для конкретного атомарного синхроимпульса цезиевого луча относительно главных синхроимпульсов Военно-морской обсерватории США. Если значение временных рядов 0, это означает, что часы цезия не потеряли и не получили времени в отношении главных часов за длительность дня. Если значение отрицательное, часы потеряли время в этот день, положительное значение означает, что часы набрали время. Для этих данных значения все отрицательные. Для некоторых приложений, таких как геодезия, важно знать, существуют ли определенные временные шкалы, где отклонение часов от ведущих часов находится на самом низком значении. Другими словами, существуют ли определенные шкалы, где часы наиболее тесно согласуются с ведущими часами? Это Преобразование Хаара полезно здесь, потому что он обладает двумя важными свойствами: Он украшает данные по шкале и разделяет энергию сигнала между шкалой.

Чтобы проиллюстрировать свойство декоррелирования, получите преобразование Хаара до уровня 6. Постройте график автокорреляционной последовательности исходных данных вместе с автокорреляцией коэффициентов вейвлета по шкале для шкал 2,4,8,16 и 32 дней. Пунктирные линии на графиках определяют 95% доверительные интервалы для входов белого шума. Значения, превышающие эти линии, указывают на значительную автокорреляцию в данных.

[s,w] = haart(clock_571,6);
helperAutoCorr(clock_571,w);

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

sigenergy = norm(clock_571,2)^2
energyByScale = cellfun(@(x)norm(x,2)^2,w);
haarenergy = norm(s,2)^2+sum(energyByScale)

sigenergy =

   2.7964e+05


haarenergy =

   2.7964e+05

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

scales = 2.^(1:6);
figure
plot(scales,energyByScale,'-o')
xlabel('Scale (days)')
set(gca,'xscale','log')
set(gca,'xtick',2.^(1:6))
ylabel('Proportion of Signal Energy')
grid on

Вы видите, что энергия минимальна для шкал 16 и 32 дня. Для вейвлета Haar (и всех вейвлетов Daubechies) вейвлет-коэффициенты в заданной шкале представляют различия между средневзвешенными средними значениями данных за длительность 1/2 длины шкалы. Этот график указывает шкалы, по которым часы цезия лучше всего согласуются с ведущими часами. Это означает, что рассмотрение данных в течение приблизительно двухнедельных или даже одномесячных периодов является более точным, чем данные в меньших или более длинных шкалах. Как упоминалось ранее, это имеет важные последствия для геодезии, где чрезвычайно точные измерения времени являются критическими.

Несмотря на то, что вейвлет Haar является прерывистым, он все еще эффективен в представлении различного рода временных рядов. Примеры включают данные подсчета и данные, где значения временных рядов связаны с некоторым определенным состоянием, которое влияет на уровень временных рядов. В качестве примера рассмотрим связь между частотой сердечных сокращений и состоянием сна.

Создайте приближения

Данные состоят из двух временных рядов. Одним из временных рядов является частота сердечных сокращений 66-дневного ребенка, отбираемого каждые 16 секунд в течение чуть более 9 часов. Временные ряды сердечных сокращений имеют целое число. Другие временные ряды являются опытным состоянием сна того же младенца за тот же период с той же частотой дискретизации. Данные о состоянии сна были оценены на основе данных EEG и EOG (движение глаз) ребенка, а не на основе частоты сердечных сокращений. Коды состояния сна: 1 = тихий сон, 2 = между тихим и активным сном, 3 = активный сон и 4 = бодрствовать. Оба временных рядов были записаны профессором Питером Флемингом, доктором Эндрю Савченко и Жанин Янг из Института здоровья детей, Королевской больницы для больных детей, Бристоль и любезно предоставлены для использования в этом примере. Постройте график данных сердечного ритма вместе с состояниями сна.

load BabyECGData
figure
yyaxis left
plot(times,HR)
ylabel('HR')
xlabel('Hrs')
YLim = [min(HR)-1 max(HR)+1];
yyaxis right
plot(times,SS)
ylabel('Sleep State')
YLim = [0.5 4.5];
title('Baby ECG and Sleep State')

Проверка данных показывает очевидную корреляцию между состоянием сна и частотой сердечных сокращений, но данные чрезвычайно зашумлены. Поскольку преобразование Хаара обеспечивает ступенчатое приближение к сигналу, это часто полезно в ситуациях, когда ответ зависит от предиктора, переменной с небольшим количеством дискретных состояний. Здесь дискретными состояниями являются четыре стадии сна. Получите приближение Haar данных пульса с помощью приближения уровня 5. Поскольку данные сердечного ритма имеют целое число, используйте 'integer' флаг для обеспечения возврата целочисленных данных. Постройте график результата.

[S,W] = haart(HR,'integer');
HaarHR = ihaart(S,W,5,'integer');
figure
hL = plotyy(times,HaarHR,times,SS);
Ax1 = hL(1);
Ax2 = hL(2);
Ax1.YLim = [min(HaarHR)-1 max(HaarHR)+1];
Ax1.YLabel.String = 'HR';
Ax2.YLim = [0.5 4.5];
Ax2.YLabel.String = 'Sleep State';
xlabel('Hrs')
title('Haar Approximation and Sleep State')

Приближение Haar более четко показывает связь между состоянием сна и данными о сердечном ритме. Вы можете оценить это изменение, посмотрев на корреляцию между необработанными данными и временными рядами состояния сна.

corr(SS,HR)

ans =

    0.5576

Теперь сравните значение 0,56 с корреляцией между данными о состоянии сна и приближением Haar

corr(SS,HaarHR)

ans =

    0.6907

Корреляция увеличилась с 0,56 до 0,69. Более совершенный вейвлет и моделирование этих данных представлены в Nason, von Sachs, & Kroisandt [1] и Nason, Sapatinas, & Sawczenko [2].

Цифровая маркировка изображений

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

Водяной знак изображение Мандрил с одним из медоносного барсука. Прочтите в изображении Mandrill. Измените его размер на 2048x2048 и отобразите результат.

coverIM = imread('mandrill.jpg');
coverIM = rgb2gray(coverIM);
coverIM = imresize(im2double(coverIM),[2048 2048]);
imagesc(coverIM)
colormap gray
title('Original Image to Watermark')
axis off
axis square

Получите преобразование Хаара изображения Мандрилла до уровня 3.

[LLorig,LHorig,HLorig,HHorig] = haart2(coverIM,3);
imagesc(LLorig)
title('Level 3 Haar Approximation')
axis off
axis square

Прочтите изображение водяного знака и измените его размер.

watermark = imread('honeybadger.jpg');
watermark = im2double(rgb2gray(watermark));
watermark = imresize(watermark,[2048 2048]);

Получите преобразование Хаара изображения водяного знака до уровня 3.

[LLw,LHw,HLw,HHw] = haart2(watermark,3);
imagesc(LLw)
colormap gray
title('Level 3 Haar Approximation--Watermark')
axis off
axis square

Добавьте водяной знак медоносного барсука к изображению Мандриля путем ослабления коэффициентов приближения уровня 3 водяного знака и вставки ослабленных коэффициентов в коэффициенты приближения уровня 3 Мэндрила.

LLwatermarked = LLorig+1e-4*LLw;
markedIM = ihaart2(LLwatermarked,LHorig,HLorig,HHorig);
imagesc(markedIM)
title('Watermarked Image')
axis off
axis square
colormap gray

Водяной знак (медоносный барсук) не отображается на изображении с водяным знаком. Поскольку вы знаете, какой алгоритм использовался для вставки водяного знака, можно восстановить водяной знак с помощью преобразования Хаара.

[LLr,LHr,HLr,HHr] = haart2(markedIM,3);
LLmarked = (LLr-LLorig).*1e4;
imagesc(LLmarked)
title('Recovered Watermark')
colormap gray
axis off
axis square