Вейвлеты и оконное преобразование Фурье для тестовых сигналов

Вейвлеты и оконное преобразование Фурье для тестовых сигналов

Упражняться и сравнивать оконное преобразование Фурье и вейвлет-спектры будем на нескольких тестовых примерах. Сгенерируйте для этих целей следующие временные ряды.

(а) Синусоида с частотой 10 Гц. Временной ряд должен иметь длину 1000 точек и частоту дискретизации 1000 Гц.

(б) Сигнал с линейно меняющейся частотой: длиной 1000 точек с частотй дискретизации 1000 Гц. Пусть ряд содержит 15 колебаний, тогда .

(в) Нестационарный сигнал, состоящий из 4-х участков: от 0 до 250мс присутствует синусоида с частотой 300 Гц, от 250мс до 500мс с частотой 200 Гц, от 500 до 750 – 100 Гц, от 750 до 1000 – 50 Гц. Частота дискретизации – те же 1000 Гц. Обеспечьте непрерывность перехода между синусоидами – следующая должна иметь начальную фазу, равную фазе на конце предыдущей синусоиды.

(г) Сигнал, представляющий собой сумму двух синусоид с частотами 500 Гц и 1 кГц. Частота дискретизации 8 кГц, длина ряда 1000 точек. К двум точкам этого ряда с номерам, например, 600 и 632 прибавьте число .

Для быстрого выполнения вейвлет-преобразования требуется использовать быстрое преобразование Фурье (БПФ). Воспользуйтесь готовой реализацией БПФ из бесплатной библиотеки Alglib (http://alglib.sources.ru/fasttransforms/fft.php). Для этого скопируйте в папку со своим проектом файлы ap.pas, fft.pas и fftbase.pas, подключите к своей программе модули ap и fft. Для работы вам потребуются следующие процедуры:

FFTR1D(Ser: TReal1DArray; N: integer; var F: TComplex1DArray) – для прямого ДПФ.

FFTR1DInv(F: TComplex1DArray; N: integer; var Ser: TReal1DArray) – для обратного ДПФ.

Здесь Ser – анализируемый временной ряд, N – количество точек в ряде, F – коэффициенты преобразования Фурье. Типы TReal1DArray и TComplex1DArray определяются как

TReal1DArray = array of double; - то есть это обычный вещественный динамический массив.

TComplex1DArray = array of Complex;

Complex = record

X, Y: double;

end;

Пример прямого и обратного преобразования для синусоиды:

Ser, Ser2: TReal1DArray;

F: TComplex1DArray;

i, N: integer;

dt: double;

begin

dt:=0.01;

N:=10000;

SetLength(Ser, N);

Ser[i]:=2 * sin(2 * Pi * i / 100);

FFTR1D(Ser, N, F);

for i:=0 to N div 2 do

Series1.AddXY(i / (N * dt), sqrt(sqr(F[i].X) + sqr(F[i].Y)));

FFTR1DInv(F, N, Ser2);

for i:=0 to N-1 do

Series2.AddXY(i * dt, Ser2[i]);

(а) Осуществите БПФ для рядов, сгенерированных вами в предыдущем задании и постройте для них спектр мощности.

(б) Создайте процедуру, осуществляющую оконное преобразование Фурье. Параметрами будут временной ряд, ширина окна и шаг, с которым оно сдвигается. Результаты записывайте в двумерный массив. По этому массиву постройте спектрограммы (зависимость амплитуд гармоник от частоты и времени (положения скользящего окна)) для всех тестовых рядов. Меняя ширину окна рассмотрите как они будет влиять на частотное и временное разрешение преобразования.

3. Вейвлет-преобразование

Рассмотрим вейвлет-преобразование (ВП) для двух материнских вейвлетов: мексиканская шляпа (MHAT):

и вейвлета Морле:

где и - параметры (традиционный выбор и ).

Их аналитически посчитанные преобразования Фурье выглядят соответственно:

и

.

(а) Постройте графики базисных вейвлетов и их спектры мощности. Проследите за тем, как эти графики будут меняться при изменении параметров a и b. По графику Фурье-образа проверьте правильность формул для перехода от масштаба вейвлет преобразования к частоте Фурье преобразования.

(б) Для вейвлета Морле рассмотрите, как будет меняться его график и Фурье-образ при изменении параметров и .

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

(г) Если максимумы на графике ВП не четкие бывает затруднительно проследить за тем, как меняется характерная частота. В этом случае строят так называемые скелетоны. Если - результат ВП, то скелетоном называют функцию вида:

Можно также использовать условия вида для получения скелетона ширины .

Если в сигнале присутствует гармоническая компонента и шум, то скелетон позволяет увидеть их раздельно. Гармонические компоненты приводят к появлению линий вытянутых по горизонтали, а шум изображается линиями, вытянутыми в основном по вертикали.

Постройте скелетоны ВП для тестовых сигналов при различном уровне шума. Например, 0, 1, 5, 10, 50, 100 и 200 процентов от стандартного отклонения сигнала.

(д) Постройте спектрограммы, ВП и их скелетоны для нескольких примеров биологических сигналов:

1) Поверхностная ЭЭГ человека (файлы eeg1.txt и eeg2.txt). Эти сигналы сняты с затылочного отведения у разных людей в состоянии т.н. пассивного бодрствования (человек расслаблен, глаза закрыты). На обеих записях просматривается -ритм – колебания с частотой порядка 8-12 Гц. Частота дискретизации сигналов 100Гц и 250Гц соответственно.