Анализ изменчивости сигналов по масштабам

Загрузить и построить график 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 дней. Для вейвлета Хаара (и всех вейвлетов Даубехия) вейвлет-коэффициенты в данной шкале представляют собой различия между средневзвешенными значениями данных в течение 1/2 длины шкалы. На этом графике показаны шкалы, по которым цезиевые часы наилучшим образом согласуются с главными часами. Это означает, что учет данных за приблизительно двухнедельные или даже одномесячные периоды является более точным, чем данные в меньших или более длинных масштабах. Как упоминалось ранее, это имеет важные последствия для геодезии, где чрезвычайно точные измерения времени имеют решающее значение.

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

Создание аппроксимаций сигналов

Данные состоят из двух временных рядов. Один временной ряд - это частота сердечных сокращений 66-дневного младенца, отбираемая каждые 16 секунд в течение чуть более 9 часов. Временной ряд частоты сердечных сокращений является целочисленным. Другой временной ряд - это опытно оцененное состояние сна одного и того же младенца за тот же период с той же частотой выборки. Данные о состоянии сна оценивались на основе данных ЭЭГ и ЭОГ ребенка (движение глаз), а не на основе частоты сердечных сокращений. Коды состояния сна: 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')

Проверка данных показывает очевидную корреляцию между состоянием сна и частотой сердечных сокращений, но данные чрезвычайно шумные. Поскольку преобразование Хаара обеспечивает приближение лестницы к сигналу, оно часто полезно в ситуациях, когда отклик зависит от прогнозирующей переменной с небольшим числом дискретных состояний. Здесь дискретными состояниями являются четыре стадии сна. Получение аппроксимации Хаара данных сердечного ритма с помощью аппроксимации уровня 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')

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

corr(SS,HR)

ans =

    0.5576

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

corr(SS,HaarHR)

ans =

    0.6907

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

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

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

Водяной знак изображения Мандрила с одним из медоносных барсуков. Читайте на изображении Мандрилла. Измените его размер на 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