Загрузите набор данных NPG2006 [1]. Данные представляют собой траекторию подповерхностного поплавка, захваченного вихрем. Постройте график смещения в восточном и северном направлениях. Треугольник обозначает исходное положение.

load npg2006
plot(npg2006.cx)
hold on
grid on
xlabel('Eastward Displacement (km)')
ylabel('Northward Displacement (km)')
plot(npg2006.cx(1),'^','markersize',11,'color','r',...
    'markerfacecolor',[1 0 0 ])

Создайте банк фильтров CWT, который можно применить к данным. Используйте вейвлет Морса по умолчанию. Период выборки данных составляет 4 часа.

fb = cwtfilterbank('SignalLength',length(npg2006.cx),'SamplingPeriod',hours(4));

Получают усредненные по времени спектры вейвлет-мощности и центральные периоды.

[tavgp,centerP] = timeSpectrum(fb,npg2006.cx);
size(tavgp)

ans = 1×3

    73     1     2

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

subplot(2,1,1)
plot(centerP,tavgp(:,1,1))
title('Counterclockwise Component')
ylabel('Power')
xlabel('Period (hrs)')
subplot(2,1,2)
plot(centerP,tavgp(:,1,2))
title('Clockwise Component')
ylabel('Power')
xlabel('Period (hrs)')

Если опустить выходные аргументы и выполнить timeSpectrum(fb,npg2006.cx) в командной строке на текущем рисунке нанесены скалограммы и усредненные по времени спектры мощности. При этом поворот поплавка по часовой стрелке фиксируется на поворотной скалограмме по часовой стрелке и усредненном по времени спектре.

Загрузить временной ряд величин солнечного магнитного поля, регистрируемых ежечасно над южным полюсом Солнца космическим аппаратом «Улисс» с 21:00 UT 4 декабря 1993 года до 12:00 UT 24 мая 1994 года. Полное описание этих данных приведено в [3] стр. 218-220. Создайте банк фильтров CWT, который можно применить к данным. Постройте график скалограммы и усредненного по времени вейвлет-спектра.

load solarMFmagnitudes
fb = cwtfilterbank('SignalLength',length(sm),'SamplingPeriod',hours(1));
timeSpectrum(fb,sm)

Получение усредненного по времени вейвлет-спектра сигнала с использованием значений по умолчанию. По умолчанию timeSpectrum нормализует мощность усредненного по времени вейвлет-спектра, чтобы равняться дисперсии сигнала. Убедитесь, что сумма спектра равна дисперсии сигнала.

tavg = timeSpectrum(fb,sm);
[var(sm) sum(tavg)]

ans = 1×2

    0.0448    0.0447

Получить усредненный по времени вейвлет-спектр сигнала, но вместо этого нормализовать мощность как функцию плотности вероятности. Убедитесь, что сумма равна 1.

tavg = timeSpectrum(fb,sm,'Normalization','pdf');
sum(tavg)

ans = 1.0000

Если установить SpectrumType кому 'density', timeSpectrum нормализует взвешенный интеграл вейвлет-спектра в соответствии со значением Normalization. Спектр имитирует функцию плотности вероятности, интеграл которой, вычисленный численно, равен значению, указанному Normalization.

Постройте график скалограммы и усредненного по времени вейвлет-спектра с типом спектра 'density' и 'pdf' нормализация.

figure
timeSpectrum(fb,sm,'SpectrumType','density','Normalization','pdf')

Для подтверждения интеграла спектра, равного 1, сначала получают усредненный по времени вейвлет-спектр с помощью 'density' тип спектра и 'pdf' нормализация.

tavg = timeSpectrum(fb,sm,'SpectrumType','density','Normalization','pdf');

По умолчанию банк фильтров использует аналитический вейвлет Морса (3,60). Получить константу допустимости для вейвлета и численно интегрировать вейвлет-спектр, используя трапециевидное правило. Помните, что CWT использует L1 нормализацию. Убедитесь, что интеграл равен 1.

ga = 3;
tbw = 60;

be = tbw/ga;
anorm = 2*exp(be/ga*(1+(log(ga)-log(be))));
cPsi = anorm^2/(2*ga).*(1/2)^(2*(be/ga)-1)*gamma(2*be/ga);

rawScales = scales(fb);
numInt = 2/cPsi*1/length(sm)*trapz(rawScales(:),tavg./rawScales(:))

numInt = 1.0000