В пакетном вейвлете выполняю декомпозицию сигнала (шумы) до шестого уровня. В последнем уровне выбираю все 64 узла. Каждый узел содержит отсчеты, возвожу каждый отсчет в квадрат и складываю. Получается энергия сигнала в данном узле. Строю график: по Х-номер узла, по Y-величина энергии.
В результате получается график зависимости энергии, приходящейся на узел, от номера узла. Частота дискретизации 8 кГц. Наблюдаю вейвлет-спектр энергии сигнала в полосе 0-4 кГц. Если разделить все 64 узла на 4 группы, то после первой и второй группы имеет место резкое падение энергии, т.е. на частотах 1 кГц и 2кГц, примерно на 3 порядка. На спектре БПФ этого не наблюдается.
Поменял вейвлет, ступени остались. Поменял шумы - остались. Взял в качестве независимого примера из Матлаба файл с шумами стиральной машины - ступени остались. Нефизично всё это.
Проверил работу кода на белом шуме (функция wgn - гаусс с равномерным распределением по частоте) - ступеней нет.
Значит программка работает правильно.
Подскажите, пожалуйста, где может быть ошибка.
PS
Рисунок можно посмотреть на
https://cloud.mail.ru/public/CX3b/Tj6XuCoZg
Вейвлет-спектр энергии шума имеет ступени
Модератор: Admin
Re: Вейвлет-спектр энергии шума имеет ступени
Частотная характеристика дерева разложения пакетного вейвлета db10.
Дерево разложения пакетного вейвлета (далее – дерево) представляет собой реализацию Коифманом и Викерхаузером усовершенствованного пирамидального алгоритма Маллата. Суть алгоритма состоит в разбиении частотного диапазона сигнала на две части: высокочастотную и низкочастотную, каждая из которых в свою очередь также разбивается на высокочастотную и низкочастотную части и так – до достижения заданного уровня.. В данном случае опустим упоминание процедуры прореживания, для данной темы она не играет роли. В результате получается дерево со сниспадающими вниз ветвями. Верхний узел, из которого начинается частотное разделение, имеет нулевой номер. Последний уровень – ряд узлов (вейвлет-фильтров) имеет номер N, означающий число делений частотных диапазонов.
Таким образом, на каждом уровне, в том числе и на последнем, ряд узлов представляет собой последовательность частотных фильтров, разделяющих весь частотный диапазон от нуля до частоты Найквиста. Это вытекает из упомянутого алгоритма Маллата и его модернизации.
Выше приведено сообщение, в котором отмечено скачкообразное изменение сигнала в узлах. Для выяснения причин этого выполнена оценка частотной характеристики (ЧХ) дерева на уровне N=6. Использованы известные функции для дерева разложения пакетного вейвлета db10 до заданного уровня. Проверка на отсутствие ошибок кода разложения выполнялась на шумовом сигнале: применил разложение с реконструкцией до заданного уровня 6 и сложил. Полученный сигнал совпал с исходным шумом, т.е. ошибок нет.
Применил указанный код для разложения синусоидального сигнала. Сначала на нулевой узел были поданы поочередно синусоиды длительностью 8с и частотами ~0.5, ~1, ~2 и ~3 кГц, с частотой дискретизации 8 кГц. Для упрощения получения оценки ЧХ в каждом из 64 узлов (фильтров на уровне N=6) отсчеты сигнала возведены в квадрат и сложены, т.е. в каждом узле (фильтре) получена энергия в полосе частот, соответствующей определенному фильтру. Полученный график можно посмотреть на https://cloud.mail.ru/public/4isX/8r6P1sqSF
Из графика следует, что на выходе дерева разложения (на уровне N=6) для каждого из перечисленных выше синусоидальных сигналов, подаваемых на нулевой узел дерева, наблюдается два сигнала, один с частотой ниже, а другой – выше частоты подаваемого сигнала.
Это противоречит пирамидальному алгоритму, согласно которому разбиению частотного диапазона сигнала производится на две неперекрывающиеся части: высокочастотную и низкочастотную.
Коллеги, как Вы думаете, чем это может быть вызвано?
Дерево разложения пакетного вейвлета (далее – дерево) представляет собой реализацию Коифманом и Викерхаузером усовершенствованного пирамидального алгоритма Маллата. Суть алгоритма состоит в разбиении частотного диапазона сигнала на две части: высокочастотную и низкочастотную, каждая из которых в свою очередь также разбивается на высокочастотную и низкочастотную части и так – до достижения заданного уровня.. В данном случае опустим упоминание процедуры прореживания, для данной темы она не играет роли. В результате получается дерево со сниспадающими вниз ветвями. Верхний узел, из которого начинается частотное разделение, имеет нулевой номер. Последний уровень – ряд узлов (вейвлет-фильтров) имеет номер N, означающий число делений частотных диапазонов.
Таким образом, на каждом уровне, в том числе и на последнем, ряд узлов представляет собой последовательность частотных фильтров, разделяющих весь частотный диапазон от нуля до частоты Найквиста. Это вытекает из упомянутого алгоритма Маллата и его модернизации.
Выше приведено сообщение, в котором отмечено скачкообразное изменение сигнала в узлах. Для выяснения причин этого выполнена оценка частотной характеристики (ЧХ) дерева на уровне N=6. Использованы известные функции для дерева разложения пакетного вейвлета db10 до заданного уровня. Проверка на отсутствие ошибок кода разложения выполнялась на шумовом сигнале: применил разложение с реконструкцией до заданного уровня 6 и сложил. Полученный сигнал совпал с исходным шумом, т.е. ошибок нет.
Применил указанный код для разложения синусоидального сигнала. Сначала на нулевой узел были поданы поочередно синусоиды длительностью 8с и частотами ~0.5, ~1, ~2 и ~3 кГц, с частотой дискретизации 8 кГц. Для упрощения получения оценки ЧХ в каждом из 64 узлов (фильтров на уровне N=6) отсчеты сигнала возведены в квадрат и сложены, т.е. в каждом узле (фильтре) получена энергия в полосе частот, соответствующей определенному фильтру. Полученный график можно посмотреть на https://cloud.mail.ru/public/4isX/8r6P1sqSF
Из графика следует, что на выходе дерева разложения (на уровне N=6) для каждого из перечисленных выше синусоидальных сигналов, подаваемых на нулевой узел дерева, наблюдается два сигнала, один с частотой ниже, а другой – выше частоты подаваемого сигнала.
Это противоречит пирамидальному алгоритму, согласно которому разбиению частотного диапазона сигнала производится на две неперекрывающиеся части: высокочастотную и низкочастотную.
Коллеги, как Вы думаете, чем это может быть вызвано?
Re: Вейвлет-спектр энергии шума имеет ступени
и частотами ~0.5, ~1, ~2 и ~3 кГц
Отрицательные частоты подразумевают, что гармонический сигнал был комплексный? Но почему частоты только отрицательные?
Применил указанный код для разложения синусоидального сигнала.
Никакого кода вы пока что не показывали.
С уважением
Александр Сергиенко
Александр Сергиенко
Re: Вейвлет-спектр энергии шума имеет ступени
Здравствуйте Александр Борисович.
1. По поводу отрицательных частот: частоту синусоиды устанавливал посредине частотного диапазона соответствующего узла, поэтому частоты получились немного отличными от 0.5, 1, 2 и 3 кГц. В связи с этим применил знак ~ указывающий на приблизительное значение частоты.
2. Попробую скомпоновать скрипт плюс файл функции и разместить здесь.
С уважением, Вениамин Дорожко.
1. По поводу отрицательных частот: частоту синусоиды устанавливал посредине частотного диапазона соответствующего узла, поэтому частоты получились немного отличными от 0.5, 1, 2 и 3 кГц. В связи с этим применил знак ~ указывающий на приблизительное значение частоты.
2. Попробую скомпоновать скрипт плюс файл функции и разместить здесь.
С уважением, Вениамин Дорожко.
Re: Вейвлет-спектр энергии шума имеет ступени
В связи с этим применил знак ~ указывающий на приблизительное значение частоты.
А, так это не минус был...
С уважением
Александр Сергиенко
Александр Сергиенко
Re: Вейвлет-спектр энергии шума имеет ступени
%% Код программы
echo on
clc
echo off
%% PARAMETERS %%%%%%%%%%%%%%%%%
Fs=8000; %% Hz - Частота дискретизации
N=2^16; %% Длина синусоиды
Ffft=0:Fs/N:Fs-Fs/N; %% Частотный вектор FFT-спектра
Nw=6; %% Заданный уровень дерева декомпозиции пакетного вейвлета
wavelet='db10'; %% Вейвлет
Fw=0:62.992:4000; %% Частотный вектор узлов на шестом уровне декомпозиции (64 узла)
%% WAVELET RESPONSE
%% ~0.5 kHz (Приблизительное значение частоты)
Fsin=0.5354*1e03; %% Точное значение частоты синусоидального сигнала
S=sin(2*pi*Fsin*(1:N)./Fs); %% Синусоидальный сигнал
[E_NODES]=my_wp_Nw6(S,Nw,wavelet); %% Функция для вычисления энергии сигнала в каждом узле
subplot(4,1,1);
plot(Ffft, abs(fft(S)), '-b', 'LineWidth',1.5); hold on;
plot(Fw, E_NODES, '-r.', 'LineWidth',1.5,'MarkerEdgeColor','r','MarkerSize',8);
axis([0 4000 0 max(abs(fft(S)))]); xlabel('F, HS');
title('fft-spectr of sinus - 0.5 kHS (blue); nodes energy of wavelet response (red)');
%% ~1.0 kHS
Fsin=1.0394*1e03;
S=sin(2*pi*Fsin*(1:N)./Fs);
[E_NODES]=my_wp_Nw6(S,Nw,wavelet);
subplot(4,1,2);
plot(Ffft, abs(fft(S)), '-b', 'LineWidth',1.5); hold on;
plot(Fw, E_NODES, '-r.', 'LineWidth',1.5,'MarkerEdgeColor','r','MarkerSize',8);
axis([0 4000 0 max(abs(fft(S)))]); xlabel('F, HS');
title('fft-spectr of sinus - 1 kHS (blue); nodes energy of wavelet response (red)');
%% ~2.0 kHS
Fsin=2.0472*1e03;
S=sin(2*pi*Fsin*(1:N)./Fs);
[E_NODES]=my_wp_Nw6(S,Nw,wavelet);
subplot(4,1,3);
plot(Ffft, abs(fft(S)), '-b', 'LineWidth',1.5); hold on;
plot(Fw, E_NODES, '-r.', 'LineWidth',1.5,'MarkerEdgeColor','r','MarkerSize',8);
axis([0 4000 0 max(abs(fft(S)))]); xlabel('F, HS');
title('fft-spectr of sinus - 2 kHS (blue); nodes energy of wavelet response (red)');
%% ~3.0 kHS
Fsin=3.0551*1e03;
S=sin(2*pi*Fsin*(1:N)./Fs);
[E_NODES]=my_wp_Nw6(S,Nw,wavelet);
subplot(4,1,4);
plot(Ffft, abs(fft(S)), '-b', 'LineWidth',1.5); hold on;
plot(Fw, E_NODES, '-r.', 'LineWidth',1.5,'MarkerEdgeColor','r','MarkerSize',8);
axis([0 4000 0 max(abs(fft(S)))]); xlabel('F, HS');
title('fft-spectr of sinus - 3 kHS (blue); nodes energy of wavelet response (red)');
%% END of SCRIPT
%% my function %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [E_NODES]=my_wp_Nw6(S, Nw,wavelet)
%% TREE Nw=6;
wpt=wpdec(S, Nw, wavelet);
%% Декомпозиция на 6-й уровень (с прореживанием) и реконструкция (восстановление до длины сигнала S)
%% Реконструкция необходима для проверки кода (см. ниже ПРОВЕРКА КОДА)
%% Для исключения реконструкции следует функцию wprcoef заменить на wpcoef , т.е. букву r исключить
%% БЕЗ РЕКОНСТРУКЦИИ РЕЗУЛЬТАТЫ ВЫЧИСЛЕНИЙ БУДУТ ТЕМИ ЖЕ
rcfs60 = wprcoef(wpt,[6 0]); rcfs610 = wprcoef(wpt,[6 10]); rcfs620 = wprcoef(wpt,[6 20]);
rcfs61 = wprcoef(wpt,[6 1]); rcfs611 = wprcoef(wpt,[6 11]); rcfs621 = wprcoef(wpt,[6 21]);
rcfs62 = wprcoef(wpt,[6 2]); rcfs612 = wprcoef(wpt,[6 12]); rcfs622 = wprcoef(wpt,[6 22]);
rcfs63 = wprcoef(wpt,[6 3]); rcfs613 = wprcoef(wpt,[6 13]); rcfs623 = wprcoef(wpt,[6 23]);
rcfs64 = wprcoef(wpt,[6 4]); rcfs614 = wprcoef(wpt,[6 14]); rcfs624 = wprcoef(wpt,[6 24]);
rcfs65 = wprcoef(wpt,[6 5]); rcfs615 = wprcoef(wpt,[6 15]); rcfs625 = wprcoef(wpt,[6 25]);
rcfs66 = wprcoef(wpt,[6 6]); rcfs616 = wprcoef(wpt,[6 16]); rcfs626 = wprcoef(wpt,[6 26]);
rcfs67 = wprcoef(wpt,[6 7]); rcfs617 = wprcoef(wpt,[6 17]); rcfs627 = wprcoef(wpt,[6 27]);
rcfs68 = wprcoef(wpt,[6 8]); rcfs618 = wprcoef(wpt,[6 18]); rcfs628 = wprcoef(wpt,[6 28]);
rcfs69 = wprcoef(wpt,[6 9]); rcfs619 = wprcoef(wpt,[6 19]); rcfs629 = wprcoef(wpt,[6 29]);
rcfs630 = wprcoef(wpt,[6 30]); rcfs640 = wprcoef(wpt,[6 40]); rcfs650 = wprcoef(wpt,[6 50]);
rcfs631 = wprcoef(wpt,[6 31]); rcfs641 = wprcoef(wpt,[6 41]); rcfs651 = wprcoef(wpt,[6 51]);
rcfs632 = wprcoef(wpt,[6 32]); rcfs642 = wprcoef(wpt,[6 42]); rcfs652 = wprcoef(wpt,[6 52]);
rcfs633 = wprcoef(wpt,[6 33]); rcfs643 = wprcoef(wpt,[6 43]); rcfs653 = wprcoef(wpt,[6 53]);
rcfs634 = wprcoef(wpt,[6 34]); rcfs644 = wprcoef(wpt,[6 44]); rcfs654 = wprcoef(wpt,[6 54]);
rcfs635 = wprcoef(wpt,[6 35]); rcfs645 = wprcoef(wpt,[6 45]); rcfs655 = wprcoef(wpt,[6 55]);
rcfs636 = wprcoef(wpt,[6 36]); rcfs646 = wprcoef(wpt,[6 46]); rcfs656 = wprcoef(wpt,[6 56]);
rcfs637 = wprcoef(wpt,[6 37]); rcfs647 = wprcoef(wpt,[6 47]); rcfs657 = wprcoef(wpt,[6 57]);
rcfs638 = wprcoef(wpt,[6 38]); rcfs648 = wprcoef(wpt,[6 48]); rcfs658 = wprcoef(wpt,[6 58]);
rcfs639 = wprcoef(wpt,[6 39]); rcfs649 = wprcoef(wpt,[6 49]); rcfs659 = wprcoef(wpt,[6 59]);
%%
rcfs660 = wprcoef(wpt,[6 60]);
rcfs661 = wprcoef(wpt,[6 61]);
rcfs662 = wprcoef(wpt,[6 62]);
rcfs663 = wprcoef(wpt,[6 63]);
%% ПРОВЕРКА КОДА
%% Суммирование реконструированного сигнала для проверки кода
% Srec=rcfs60 +rcfs61 +rcfs62 +rcfs63 +rcfs64 +rcfs65 +rcfs66 +rcfs67 +rcfs68 +rcfs69+...
% rcfs610 +rcfs611 +rcfs612 +rcfs613 +rcfs614 +rcfs615 +rcfs616 +rcfs617 +rcfs618 +rcfs619+...
% rcfs620 +rcfs621 +rcfs622 +rcfs623 +rcfs624 +rcfs625 +rcfs626 +rcfs627 +rcfs628 +rcfs629+...
% rcfs630 +rcfs631 +rcfs632 +rcfs633 +rcfs634 +rcfs635 +rcfs636 +rcfs637 +rcfs638 +rcfs639+...
% rcfs640 +rcfs641 +rcfs642 +rcfs643 +rcfs644 +rcfs645 +rcfs646 +rcfs647 +rcfs648 +rcfs649+...
% rcfs650 +rcfs651 +rcfs652 +rcfs653 +rcfs654 +rcfs655 +rcfs656 +rcfs657 +rcfs658 +rcfs659+...
% rcfs660 +rcfs661 +rcfs662 +rcfs663 ;
%
% Delta=S-Srec; plot(delta);
%
%% Сравнение реконструированного сигнала с исходной синусоидой.
%% Разность Delta на графике составляет примерно 1e-011 ,
%% что вполне достаточно для вейвлет декомпозиции плюс реконструкции
%% т.е. КОД РАБОТАЕТ ВЕРНО
%% ПРИМЕЧАНИЕ: без реконструкции сигналов (вейвлет-коэффициентов в узлах)
%% вычисляемые графики будут теми же
%% Энергия сигналов в каждом узле
E_NODES=[sum( rcfs60.^2) sum( rcfs61.^2) sum( rcfs62.^2) sum( rcfs63.^2) sum( rcfs64.^2)...
sum( rcfs65.^2) sum( rcfs66.^2) sum( rcfs67.^2) sum( rcfs68.^2) sum( rcfs69.^2)...
sum( rcfs610.^2) sum( rcfs611.^2) sum( rcfs612.^2) sum( rcfs613.^2) sum( rcfs614.^2)...
sum( rcfs615.^2) sum( rcfs616.^2) sum( rcfs617.^2) sum( rcfs618.^2) sum( rcfs619.^2)...
sum( rcfs620.^2) sum( rcfs621.^2) sum( rcfs622.^2) sum( rcfs623.^2) sum( rcfs624.^2)...
sum( rcfs625.^2) sum( rcfs626.^2) sum( rcfs627.^2) sum( rcfs628.^2) sum( rcfs629.^2)...
sum( rcfs630.^2) sum( rcfs631.^2) sum( rcfs632.^2) sum( rcfs633.^2) sum( rcfs634.^2)...
sum( rcfs635.^2) sum( rcfs636.^2) sum( rcfs637.^2) sum( rcfs638.^2) sum( rcfs639.^2)...
sum( rcfs640.^2) sum( rcfs641.^2) sum( rcfs642.^2) sum( rcfs643.^2) sum( rcfs644.^2)...
sum( rcfs645.^2) sum( rcfs646.^2) sum( rcfs647.^2) sum( rcfs648.^2) sum( rcfs649.^2)...
sum( rcfs650.^2) sum( rcfs651.^2) sum( rcfs652.^2) sum( rcfs653.^2) sum( rcfs654.^2)...
sum( rcfs655.^2) sum( rcfs656.^2) sum( rcfs657.^2) sum( rcfs658.^2) sum( rcfs659.^2)...
sum( rcfs660.^2) sum( rcfs661.^2) sum( rcfs662.^2) sum( rcfs663.^2) ];
%% END
echo on
clc
echo off
%% PARAMETERS %%%%%%%%%%%%%%%%%
Fs=8000; %% Hz - Частота дискретизации
N=2^16; %% Длина синусоиды
Ffft=0:Fs/N:Fs-Fs/N; %% Частотный вектор FFT-спектра
Nw=6; %% Заданный уровень дерева декомпозиции пакетного вейвлета
wavelet='db10'; %% Вейвлет
Fw=0:62.992:4000; %% Частотный вектор узлов на шестом уровне декомпозиции (64 узла)
%% WAVELET RESPONSE
%% ~0.5 kHz (Приблизительное значение частоты)
Fsin=0.5354*1e03; %% Точное значение частоты синусоидального сигнала
S=sin(2*pi*Fsin*(1:N)./Fs); %% Синусоидальный сигнал
[E_NODES]=my_wp_Nw6(S,Nw,wavelet); %% Функция для вычисления энергии сигнала в каждом узле
subplot(4,1,1);
plot(Ffft, abs(fft(S)), '-b', 'LineWidth',1.5); hold on;
plot(Fw, E_NODES, '-r.', 'LineWidth',1.5,'MarkerEdgeColor','r','MarkerSize',8);
axis([0 4000 0 max(abs(fft(S)))]); xlabel('F, HS');
title('fft-spectr of sinus - 0.5 kHS (blue); nodes energy of wavelet response (red)');
%% ~1.0 kHS
Fsin=1.0394*1e03;
S=sin(2*pi*Fsin*(1:N)./Fs);
[E_NODES]=my_wp_Nw6(S,Nw,wavelet);
subplot(4,1,2);
plot(Ffft, abs(fft(S)), '-b', 'LineWidth',1.5); hold on;
plot(Fw, E_NODES, '-r.', 'LineWidth',1.5,'MarkerEdgeColor','r','MarkerSize',8);
axis([0 4000 0 max(abs(fft(S)))]); xlabel('F, HS');
title('fft-spectr of sinus - 1 kHS (blue); nodes energy of wavelet response (red)');
%% ~2.0 kHS
Fsin=2.0472*1e03;
S=sin(2*pi*Fsin*(1:N)./Fs);
[E_NODES]=my_wp_Nw6(S,Nw,wavelet);
subplot(4,1,3);
plot(Ffft, abs(fft(S)), '-b', 'LineWidth',1.5); hold on;
plot(Fw, E_NODES, '-r.', 'LineWidth',1.5,'MarkerEdgeColor','r','MarkerSize',8);
axis([0 4000 0 max(abs(fft(S)))]); xlabel('F, HS');
title('fft-spectr of sinus - 2 kHS (blue); nodes energy of wavelet response (red)');
%% ~3.0 kHS
Fsin=3.0551*1e03;
S=sin(2*pi*Fsin*(1:N)./Fs);
[E_NODES]=my_wp_Nw6(S,Nw,wavelet);
subplot(4,1,4);
plot(Ffft, abs(fft(S)), '-b', 'LineWidth',1.5); hold on;
plot(Fw, E_NODES, '-r.', 'LineWidth',1.5,'MarkerEdgeColor','r','MarkerSize',8);
axis([0 4000 0 max(abs(fft(S)))]); xlabel('F, HS');
title('fft-spectr of sinus - 3 kHS (blue); nodes energy of wavelet response (red)');
%% END of SCRIPT
%% my function %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [E_NODES]=my_wp_Nw6(S, Nw,wavelet)
%% TREE Nw=6;
wpt=wpdec(S, Nw, wavelet);
%% Декомпозиция на 6-й уровень (с прореживанием) и реконструкция (восстановление до длины сигнала S)
%% Реконструкция необходима для проверки кода (см. ниже ПРОВЕРКА КОДА)
%% Для исключения реконструкции следует функцию wprcoef заменить на wpcoef , т.е. букву r исключить
%% БЕЗ РЕКОНСТРУКЦИИ РЕЗУЛЬТАТЫ ВЫЧИСЛЕНИЙ БУДУТ ТЕМИ ЖЕ
rcfs60 = wprcoef(wpt,[6 0]); rcfs610 = wprcoef(wpt,[6 10]); rcfs620 = wprcoef(wpt,[6 20]);
rcfs61 = wprcoef(wpt,[6 1]); rcfs611 = wprcoef(wpt,[6 11]); rcfs621 = wprcoef(wpt,[6 21]);
rcfs62 = wprcoef(wpt,[6 2]); rcfs612 = wprcoef(wpt,[6 12]); rcfs622 = wprcoef(wpt,[6 22]);
rcfs63 = wprcoef(wpt,[6 3]); rcfs613 = wprcoef(wpt,[6 13]); rcfs623 = wprcoef(wpt,[6 23]);
rcfs64 = wprcoef(wpt,[6 4]); rcfs614 = wprcoef(wpt,[6 14]); rcfs624 = wprcoef(wpt,[6 24]);
rcfs65 = wprcoef(wpt,[6 5]); rcfs615 = wprcoef(wpt,[6 15]); rcfs625 = wprcoef(wpt,[6 25]);
rcfs66 = wprcoef(wpt,[6 6]); rcfs616 = wprcoef(wpt,[6 16]); rcfs626 = wprcoef(wpt,[6 26]);
rcfs67 = wprcoef(wpt,[6 7]); rcfs617 = wprcoef(wpt,[6 17]); rcfs627 = wprcoef(wpt,[6 27]);
rcfs68 = wprcoef(wpt,[6 8]); rcfs618 = wprcoef(wpt,[6 18]); rcfs628 = wprcoef(wpt,[6 28]);
rcfs69 = wprcoef(wpt,[6 9]); rcfs619 = wprcoef(wpt,[6 19]); rcfs629 = wprcoef(wpt,[6 29]);
rcfs630 = wprcoef(wpt,[6 30]); rcfs640 = wprcoef(wpt,[6 40]); rcfs650 = wprcoef(wpt,[6 50]);
rcfs631 = wprcoef(wpt,[6 31]); rcfs641 = wprcoef(wpt,[6 41]); rcfs651 = wprcoef(wpt,[6 51]);
rcfs632 = wprcoef(wpt,[6 32]); rcfs642 = wprcoef(wpt,[6 42]); rcfs652 = wprcoef(wpt,[6 52]);
rcfs633 = wprcoef(wpt,[6 33]); rcfs643 = wprcoef(wpt,[6 43]); rcfs653 = wprcoef(wpt,[6 53]);
rcfs634 = wprcoef(wpt,[6 34]); rcfs644 = wprcoef(wpt,[6 44]); rcfs654 = wprcoef(wpt,[6 54]);
rcfs635 = wprcoef(wpt,[6 35]); rcfs645 = wprcoef(wpt,[6 45]); rcfs655 = wprcoef(wpt,[6 55]);
rcfs636 = wprcoef(wpt,[6 36]); rcfs646 = wprcoef(wpt,[6 46]); rcfs656 = wprcoef(wpt,[6 56]);
rcfs637 = wprcoef(wpt,[6 37]); rcfs647 = wprcoef(wpt,[6 47]); rcfs657 = wprcoef(wpt,[6 57]);
rcfs638 = wprcoef(wpt,[6 38]); rcfs648 = wprcoef(wpt,[6 48]); rcfs658 = wprcoef(wpt,[6 58]);
rcfs639 = wprcoef(wpt,[6 39]); rcfs649 = wprcoef(wpt,[6 49]); rcfs659 = wprcoef(wpt,[6 59]);
%%
rcfs660 = wprcoef(wpt,[6 60]);
rcfs661 = wprcoef(wpt,[6 61]);
rcfs662 = wprcoef(wpt,[6 62]);
rcfs663 = wprcoef(wpt,[6 63]);
%% ПРОВЕРКА КОДА
%% Суммирование реконструированного сигнала для проверки кода
% Srec=rcfs60 +rcfs61 +rcfs62 +rcfs63 +rcfs64 +rcfs65 +rcfs66 +rcfs67 +rcfs68 +rcfs69+...
% rcfs610 +rcfs611 +rcfs612 +rcfs613 +rcfs614 +rcfs615 +rcfs616 +rcfs617 +rcfs618 +rcfs619+...
% rcfs620 +rcfs621 +rcfs622 +rcfs623 +rcfs624 +rcfs625 +rcfs626 +rcfs627 +rcfs628 +rcfs629+...
% rcfs630 +rcfs631 +rcfs632 +rcfs633 +rcfs634 +rcfs635 +rcfs636 +rcfs637 +rcfs638 +rcfs639+...
% rcfs640 +rcfs641 +rcfs642 +rcfs643 +rcfs644 +rcfs645 +rcfs646 +rcfs647 +rcfs648 +rcfs649+...
% rcfs650 +rcfs651 +rcfs652 +rcfs653 +rcfs654 +rcfs655 +rcfs656 +rcfs657 +rcfs658 +rcfs659+...
% rcfs660 +rcfs661 +rcfs662 +rcfs663 ;
%
% Delta=S-Srec; plot(delta);
%
%% Сравнение реконструированного сигнала с исходной синусоидой.
%% Разность Delta на графике составляет примерно 1e-011 ,
%% что вполне достаточно для вейвлет декомпозиции плюс реконструкции
%% т.е. КОД РАБОТАЕТ ВЕРНО
%% ПРИМЕЧАНИЕ: без реконструкции сигналов (вейвлет-коэффициентов в узлах)
%% вычисляемые графики будут теми же
%% Энергия сигналов в каждом узле
E_NODES=[sum( rcfs60.^2) sum( rcfs61.^2) sum( rcfs62.^2) sum( rcfs63.^2) sum( rcfs64.^2)...
sum( rcfs65.^2) sum( rcfs66.^2) sum( rcfs67.^2) sum( rcfs68.^2) sum( rcfs69.^2)...
sum( rcfs610.^2) sum( rcfs611.^2) sum( rcfs612.^2) sum( rcfs613.^2) sum( rcfs614.^2)...
sum( rcfs615.^2) sum( rcfs616.^2) sum( rcfs617.^2) sum( rcfs618.^2) sum( rcfs619.^2)...
sum( rcfs620.^2) sum( rcfs621.^2) sum( rcfs622.^2) sum( rcfs623.^2) sum( rcfs624.^2)...
sum( rcfs625.^2) sum( rcfs626.^2) sum( rcfs627.^2) sum( rcfs628.^2) sum( rcfs629.^2)...
sum( rcfs630.^2) sum( rcfs631.^2) sum( rcfs632.^2) sum( rcfs633.^2) sum( rcfs634.^2)...
sum( rcfs635.^2) sum( rcfs636.^2) sum( rcfs637.^2) sum( rcfs638.^2) sum( rcfs639.^2)...
sum( rcfs640.^2) sum( rcfs641.^2) sum( rcfs642.^2) sum( rcfs643.^2) sum( rcfs644.^2)...
sum( rcfs645.^2) sum( rcfs646.^2) sum( rcfs647.^2) sum( rcfs648.^2) sum( rcfs649.^2)...
sum( rcfs650.^2) sum( rcfs651.^2) sum( rcfs652.^2) sum( rcfs653.^2) sum( rcfs654.^2)...
sum( rcfs655.^2) sum( rcfs656.^2) sum( rcfs657.^2) sum( rcfs658.^2) sum( rcfs659.^2)...
sum( rcfs660.^2) sum( rcfs661.^2) sum( rcfs662.^2) sum( rcfs663.^2) ];
%% END
Re: Вейвлет-спектр энергии шума имеет ступени
К сообщению (см. выше) о частотной характеристике дерева разложения пакетного вейвлета или почему не получается частотная характеристика пакетной вейвлет-декомпозиции.
Рассмотрим вейвлет-декомпозицию пакетным вейвлетом 'db40' сигнала Z (белый шум) длительностью 5 секунд с частотой дискретизации 10000 Гц (см. Приведенный ниже код). Результаты декомпозиции на трех уровнях представлены на рисунке (рисунок можно посмотреть на https://cloud.mail.ru/public/EduB/T4huqCgGi )
Известно, что в дереве пакетной вейвлет-декомпозиции (fig=drawtree(T)) узлы пронумерованы следующим образом: [Nw n], где Nw - уровень декомпозиции, n - порядковый номер узла, исчисляемый слева-направо на заданном уровне.
Согласно алгоритму пирамидальной декомпозиции каждый узел (на заданном уровне) занимает свой частотный диапазон. С увеличением порядкового номера (n), на заданном уровне, частотный диапазон узла сдвигается по оси частоты в сторону увеличения частот. Рассматривая несколько узлов в порядке возрастания их номеров, мы ожидаем, что fft-спектры их вейвлет-коэффициентов (частотные диапазоны) также будут располагаться последовательно, один за другим, по оси частот.
Например, на 1-ом уровне (см. рисунок - адрес см. выше) узел [1 0] соответствует диапазону 0-2500 Гц, а узел [1 1] – диапазону 2500-5000 Гц. Здесь 5000 Гц это частота Найквиста при частоте дискретизации 10000 Гц. Т.е. сначала идет частотный диапазон узла [1 0], потом – следующий за ним частотный диапазон узла [1 1].
Из рисунка следует, что на 1-ом уровне соблюдается указанный порядок.
На 2-ом уровне (см. рисунок) декомпозиции указанный порядок нарушается. Узлы [2 2] и [2 3] в дереве декомпозиции идут один за другим, а fft-спектры их вейвлет-коэффициентов имеют обратный порядок: сначала по частоте располагается спектр узла [2 3], а затем – узла [2 2].
Аналогичное нарушение порядка расположения частотных диапазонов наблюдается на 3-ем уровне.
Сначала по оси частот наблюдается fft-спектр узла [3 3], за ним – узла [3 2]. Аналогичное реверсное расположение частотных диапазонов имеет место и для узлов [3 5] и [3 4].
Кроме того, имеет место сдвиг частот узлов [3 4] и [3 5] в сторону увеличения, который привел к образованию пустого частотного диапазона 3000-3500 Гц.
Примечание: в заголовках к рисункам каждого уровня декомпозиции указаны номера узлов и соответствующие им цвета графиков.
Из вышесказанного следует, невозможность получить АЧХ по-видимому обусловлена следующими особенностями пакетной вейвлет-декомпозиции:
- расположение частотных диапазонов может не соответствовать порядковому номеру их узлов;
- частотные диапазоны узлов могут не соответствовать их расчетным значениям (сдвигаются);
- имеет место образование пустых частотных диапазонов;
- функция дерева пакетной вейвлет-декомпозиции T=wpdec(Z, Nw, ‘wavelet’, ‘entropy’) имеет параметр энтропии (‘entropy’), роль которого в дереве декомпозиции в литературе не раскрывается.
%% Ниже приводится код программы вейвлет-декомпозиции сигнала «белый шум».
wavelet='db40';
Fs=10000; %% Hz
Z=1*wgn(1, 5*Fs,-20);%% white noise
N=length(Z);
Ffft=0:Fs/N:Fs-Fs/N; %% frequency, Hz
time=(0:1:N-1)/Fs; %% time, s
Y00=abs(fft(Z)/N); %% FFT-specter of white noise
Nw=3; %% Max Level of wavelet-decomposition
T=wpdec(Z, Nw, wavelet, 'shannon'); %% fig=drawtree(T);
subplot(4,1,1);
plot(Ffft, Y00, '-k'); axis([0 Fs/2, 0 max(Y00)]);
title('FFT-specter of white noise');
%% LEVEL=1
rcfs10 = wprcoef(T,[1 0]); Y10=abs(fft(rcfs10)/N); %% reconstructed wavelet-coefficients
rcfs11 = wprcoef(T,[1 1]); Y11=abs(fft(rcfs11)/N);
subplot(4,1,2);
plot(Ffft, Y10, '-k'); hold on;
title('FFT-specter of wavelet-coefficients on Level=1; NodeColor: [1 0]-k, [1 1]-b');
plot(Ffft, Y11, '-b'); axis([0 Fs/2, 0 max(Y10)]);
%% LEVEL=2
rcfs20 = wprcoef(T,[2 0]); Y20=abs(fft(rcfs20)/N); %% reconstructed wavelet-coefficients
rcfs21 = wprcoef(T,[2 1]); Y21=abs(fft(rcfs21)/N);
rcfs22 = wprcoef(T,[2 2]); Y22=abs(fft(rcfs22)/N);
rcfs23 = wprcoef(T,[2 3]); Y23=abs(fft(rcfs23)/N);
subplot(4,1,3);
plot(Ffft, Y20, '-k'); hold on;
title('FFT-specter of wavelet-coefficients on Level=2; NodeColor: [2 0]-k, [2 1]-b, [2 2]-r, [2 3]-m');
plot(Ffft, Y21, '-b');
plot(Ffft, Y22, '-r');
plot(Ffft, Y23, '-m');
axis([0 Fs/2, 0 max(Y20)]);
%% LEVEL=3
rcfs30 = wprcoef(T,[3 0]); Y30=abs(fft(rcfs30)/N); %% reconstructed wavelet-coefficients
rcfs31 = wprcoef(T,[3 1]); Y31=abs(fft(rcfs31)/N);
rcfs32 = wprcoef(T,[3 2]); Y32=abs(fft(rcfs32)/N);
rcfs33 = wprcoef(T,[3 3]); Y33=abs(fft(rcfs33)/N);
rcfs34 = wprcoef(T,[3 4]); Y34=abs(fft(rcfs34)/N);
rcfs35 = wprcoef(T,[3 5]); Y35=abs(fft(rcfs35)/N);
subplot(4,1,4);
plot(Ffft, Y30, '-k'); hold on;
title('FFT-specter of wavelet-coefficients on Level=3; NodeColor: [3 0]-k, [3 1]-b, [3 2]-r, [3 3]-m, [3 4]-g, [3 5]-c');
plot(Ffft, Y31, '-b');
plot(Ffft, Y32, '-r');
plot(Ffft, Y33, '-m');
plot(Ffft, Y34, '-g');
plot(Ffft, Y35, '-c');
axis([0 Fs/2, 0 max(Y30)]); xlabel('F, Hz');
%% END OF CODE
Рассмотрим вейвлет-декомпозицию пакетным вейвлетом 'db40' сигнала Z (белый шум) длительностью 5 секунд с частотой дискретизации 10000 Гц (см. Приведенный ниже код). Результаты декомпозиции на трех уровнях представлены на рисунке (рисунок можно посмотреть на https://cloud.mail.ru/public/EduB/T4huqCgGi )
Известно, что в дереве пакетной вейвлет-декомпозиции (fig=drawtree(T)) узлы пронумерованы следующим образом: [Nw n], где Nw - уровень декомпозиции, n - порядковый номер узла, исчисляемый слева-направо на заданном уровне.
Согласно алгоритму пирамидальной декомпозиции каждый узел (на заданном уровне) занимает свой частотный диапазон. С увеличением порядкового номера (n), на заданном уровне, частотный диапазон узла сдвигается по оси частоты в сторону увеличения частот. Рассматривая несколько узлов в порядке возрастания их номеров, мы ожидаем, что fft-спектры их вейвлет-коэффициентов (частотные диапазоны) также будут располагаться последовательно, один за другим, по оси частот.
Например, на 1-ом уровне (см. рисунок - адрес см. выше) узел [1 0] соответствует диапазону 0-2500 Гц, а узел [1 1] – диапазону 2500-5000 Гц. Здесь 5000 Гц это частота Найквиста при частоте дискретизации 10000 Гц. Т.е. сначала идет частотный диапазон узла [1 0], потом – следующий за ним частотный диапазон узла [1 1].
Из рисунка следует, что на 1-ом уровне соблюдается указанный порядок.
На 2-ом уровне (см. рисунок) декомпозиции указанный порядок нарушается. Узлы [2 2] и [2 3] в дереве декомпозиции идут один за другим, а fft-спектры их вейвлет-коэффициентов имеют обратный порядок: сначала по частоте располагается спектр узла [2 3], а затем – узла [2 2].
Аналогичное нарушение порядка расположения частотных диапазонов наблюдается на 3-ем уровне.
Сначала по оси частот наблюдается fft-спектр узла [3 3], за ним – узла [3 2]. Аналогичное реверсное расположение частотных диапазонов имеет место и для узлов [3 5] и [3 4].
Кроме того, имеет место сдвиг частот узлов [3 4] и [3 5] в сторону увеличения, который привел к образованию пустого частотного диапазона 3000-3500 Гц.
Примечание: в заголовках к рисункам каждого уровня декомпозиции указаны номера узлов и соответствующие им цвета графиков.
Из вышесказанного следует, невозможность получить АЧХ по-видимому обусловлена следующими особенностями пакетной вейвлет-декомпозиции:
- расположение частотных диапазонов может не соответствовать порядковому номеру их узлов;
- частотные диапазоны узлов могут не соответствовать их расчетным значениям (сдвигаются);
- имеет место образование пустых частотных диапазонов;
- функция дерева пакетной вейвлет-декомпозиции T=wpdec(Z, Nw, ‘wavelet’, ‘entropy’) имеет параметр энтропии (‘entropy’), роль которого в дереве декомпозиции в литературе не раскрывается.
%% Ниже приводится код программы вейвлет-декомпозиции сигнала «белый шум».
wavelet='db40';
Fs=10000; %% Hz
Z=1*wgn(1, 5*Fs,-20);%% white noise
N=length(Z);
Ffft=0:Fs/N:Fs-Fs/N; %% frequency, Hz
time=(0:1:N-1)/Fs; %% time, s
Y00=abs(fft(Z)/N); %% FFT-specter of white noise
Nw=3; %% Max Level of wavelet-decomposition
T=wpdec(Z, Nw, wavelet, 'shannon'); %% fig=drawtree(T);
subplot(4,1,1);
plot(Ffft, Y00, '-k'); axis([0 Fs/2, 0 max(Y00)]);
title('FFT-specter of white noise');
%% LEVEL=1
rcfs10 = wprcoef(T,[1 0]); Y10=abs(fft(rcfs10)/N); %% reconstructed wavelet-coefficients
rcfs11 = wprcoef(T,[1 1]); Y11=abs(fft(rcfs11)/N);
subplot(4,1,2);
plot(Ffft, Y10, '-k'); hold on;
title('FFT-specter of wavelet-coefficients on Level=1; NodeColor: [1 0]-k, [1 1]-b');
plot(Ffft, Y11, '-b'); axis([0 Fs/2, 0 max(Y10)]);
%% LEVEL=2
rcfs20 = wprcoef(T,[2 0]); Y20=abs(fft(rcfs20)/N); %% reconstructed wavelet-coefficients
rcfs21 = wprcoef(T,[2 1]); Y21=abs(fft(rcfs21)/N);
rcfs22 = wprcoef(T,[2 2]); Y22=abs(fft(rcfs22)/N);
rcfs23 = wprcoef(T,[2 3]); Y23=abs(fft(rcfs23)/N);
subplot(4,1,3);
plot(Ffft, Y20, '-k'); hold on;
title('FFT-specter of wavelet-coefficients on Level=2; NodeColor: [2 0]-k, [2 1]-b, [2 2]-r, [2 3]-m');
plot(Ffft, Y21, '-b');
plot(Ffft, Y22, '-r');
plot(Ffft, Y23, '-m');
axis([0 Fs/2, 0 max(Y20)]);
%% LEVEL=3
rcfs30 = wprcoef(T,[3 0]); Y30=abs(fft(rcfs30)/N); %% reconstructed wavelet-coefficients
rcfs31 = wprcoef(T,[3 1]); Y31=abs(fft(rcfs31)/N);
rcfs32 = wprcoef(T,[3 2]); Y32=abs(fft(rcfs32)/N);
rcfs33 = wprcoef(T,[3 3]); Y33=abs(fft(rcfs33)/N);
rcfs34 = wprcoef(T,[3 4]); Y34=abs(fft(rcfs34)/N);
rcfs35 = wprcoef(T,[3 5]); Y35=abs(fft(rcfs35)/N);
subplot(4,1,4);
plot(Ffft, Y30, '-k'); hold on;
title('FFT-specter of wavelet-coefficients on Level=3; NodeColor: [3 0]-k, [3 1]-b, [3 2]-r, [3 3]-m, [3 4]-g, [3 5]-c');
plot(Ffft, Y31, '-b');
plot(Ffft, Y32, '-r');
plot(Ffft, Y33, '-m');
plot(Ffft, Y34, '-g');
plot(Ffft, Y35, '-c');
axis([0 Fs/2, 0 max(Y30)]); xlabel('F, Hz');
%% END OF CODE
Re: Вейвлет-спектр энергии шума имеет ступени
Benjamin писал(а):- имеет место образование пустых частотных диапазонов;
.....
rcfs34 = wprcoef(T,[3 4]); Y34=abs(fft(rcfs34)/N);
rcfs35 = wprcoef(T,[3 5]); Y35=abs(fft(rcfs35)/N);
Имеет место невнимательность. Пропущены еще два диапазона (3-й уровень - 2^3=8 диапазонов)
rcfs36 = wprcoef(T,[3 6]); Y36=abs(fft(rcfs36)/N);
rcfs37 = wprcoef(T,[3 7]); Y37=abs(fft(rcfs37)/N);
Re: Вейвлет-спектр энергии шума имеет ступени
Добрый день, коллеги.
Прежде всего, хочу поблагодарить коллегу U235 за внимание к вопросу и высказанное замечание.
Мне не удалось до сих пор выяснить, почему порядок следования по частоте спектров вейвлет-коэффициентов не соответствует порядку возрастания номеров их узлов последнего уровня в дереве декомпозиции сигнала пакетным вейвлетом.
Ниже предлагаю Вашему вниманию изложение сути вопроса.
Дерево декомпозиции сигнала пакетным вейвлетом представляет собой реализацию Кйфманом и Викерхаузером усовершенствованного пирамидального алгоритма Маллата. Суть алгоритма состоит в разбиении частотного диапазона сигнала на две части: высокочастотную и низкочастотную, каждая из которых в свою очередь также разбивается на высокочастотную и низкочастотную части и так – до достижения заданного уровня. В результате получается дерево со сниспадающими вниз ветвями. Верхний узел, из которого начинается частотное разделение, имеет нулевой номер. Последний уровень имеет номер Nw, означающий число делений частотных диапазонов. Таким образом, на каждом уровне, в том числе и на последнем, ряд узлов можно представить как последовательность вейвлет-спектров, спектры которых разделяют весь частотный диапазон от нуля до частоты Найквиста.
Узлы обозначаются следующим образом [Nw, n], где Nw – уровень декомпозиции, n – порядковый номер узла на заданном уровне Nw.
Ожидается, что спектры вейвлет-коэффициентов узлов идут по частоте один за другим по мере увеличения порядкового номера узла.
В данном случае прореженные вейвлет-коэффициенты подвергаются процедуре реконструкции для восстановления числа отсчетов. Это делается для того, чтобы спектры вейвлет-коэффициентов соответствовали частотному диапазону исходного сигнала. Тогда спектр вейвлет-коэффициентов сигнала начального узла [Nw, 0] будет находиться в начале частотного диапазона, а спектр вейвлет-коэффициентов сигнала конечного узла [Nw, 2^Nw-1] будет находиться в конце частотного диапазона, примыкать к частоте 5 кГц (в данном примере частота дискретизации Fs равна 10 кГц).
Ниже приводится код, который позволяет получить диаграмму спектров в осях: ордината – номера узлов, абсцисса – частота спектров. Уровень декомпозиции Nw=5.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
echo on
clc
echo off
%% Сигнал
Fs=10000; % Частота дискретизации
S=wgn(1, 50*Fs, 0); % Белый шум
N=length(S); % Длина вектора сигнала
Ffft=0:Fs/N:Fs-Fs/N; % Частота для спектра Фурье
%% WP
Nw=5; %% Уровень декомпозиции
wavelet = 'db40'; % Вейвлет
T = wptree(2,Nw,S,wavelet); % Дерево декомпозиции пакетного вейвлета
%% Цикл для вычисления схемы частотного положения спектров
% вейвлет-коэффициентов на последнем уровне Nw=5
for i=0:1:2^Nw-1 %% i - порядковый номер узла на заданном уровне Nw
rcfs = wprcoef(T,[Nw i]); % Вейвлет-коэффициенты в заданном узле i
Yrcfs=abs(fft(rcfs)); % Спектр фурье вейвлет-коэффициентов
% Так как спектры оказались сгруппированными по четыре, то решено
% выполнить их окраску «по четверкам»
% Цвета спектров каждой четверки:
% ЧЕРНЫЙ – узел с номером [5, 0], т.е. с номером 0,
% а также каждый следующий узел с номером «плюс 4»
% ГОЛУБОЙ – узел с номером 1 и далее – каждый с номером «плюс 4»
% ЗЕЛЕНЫЙ – узел с номером 2 и далее – каждый с номером «плюс 4»
% КРАСНЫЙ – узел с номером 3 и далее – каждый с номером «плюс 4»
% Цветовая дифференциация позволяет обратить внимание на инверсию
% порядка следования третьего (зеленый) и четвертого (красный)
% спектра внутри каждой четверки.
% И не только, имеет место инверсия всей четверки.
% Зачем, не понятно
vect32=['k' 'b' 'g' 'r' ...
'k' 'b' 'g' 'r' ...
'k' 'b' 'g' 'r' ...
'k' 'b' 'g' 'r' ...
'k' 'b' 'g' 'r' ...
'k' 'b' 'g' 'r' ...
'k' 'b' 'g' 'r' ...
'k' 'b' 'g' 'r'];
vect=[vect32 vect32 vect32 vect32];
%% PLOT
plot(Ffft, i+Yrcfs./max(Yrcfs), vect(i+1)); grid on; hold on; title('Nw=5; wavelet=db40');
xlabel('F, Hz'); ylabel('Node number');
set(gca,'FontName','TimesET'); set(gca,'FontWeight', 'normal');set(gca,'FontSize',14);
axis([0 5000 0 i+1]);
end
%% END %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Рисунок, полученный с помощью этого кода, размещен по адресу:
http://www.imageup.ru/img2/3408073/wp_w ... y.jpg.html
На рисунке над спектрами для улучшения чтения дополнительно проставлены номера узлов, которые соответствуют номерам на оси ординат.
Из рисунка видно, что спектры узлов не располагаются по оси частот последовательно, один за другим, сообразно номерам их узлов. Самым крайним спектром (справа, вблизи частоты 5 кГц) оказался спектр узла с номером 16, а спектр последнего узла (с номером 31) оказался в средине частотного диапазона – 3300-3500 Гц.
Если исходить из трактовки усовершенствованного пирамидального алгоритма Мала (Суть алгоритма состоит в разбиении частотного диапазона сигнала на две части: высокочастотную и низкочастотную, каждая из которых в свою очередь также разбивается на высокочастотную и низкочастотную части и так – до достижения заданного уровня), то частотные диапазоны спектров узлов должны идти один за другим по мере увеличения порядкового номера узла.
Однако, все перепутано местами.
Коллеги, подскажите, пожалуйста, как добиться, чтобы спектры вейвлет-коэффициентов узлов шли по частоте один за другим соответственно номерам их узлов.
С уважением, Вениамин.
Прежде всего, хочу поблагодарить коллегу U235 за внимание к вопросу и высказанное замечание.
Мне не удалось до сих пор выяснить, почему порядок следования по частоте спектров вейвлет-коэффициентов не соответствует порядку возрастания номеров их узлов последнего уровня в дереве декомпозиции сигнала пакетным вейвлетом.
Ниже предлагаю Вашему вниманию изложение сути вопроса.
Дерево декомпозиции сигнала пакетным вейвлетом представляет собой реализацию Кйфманом и Викерхаузером усовершенствованного пирамидального алгоритма Маллата. Суть алгоритма состоит в разбиении частотного диапазона сигнала на две части: высокочастотную и низкочастотную, каждая из которых в свою очередь также разбивается на высокочастотную и низкочастотную части и так – до достижения заданного уровня. В результате получается дерево со сниспадающими вниз ветвями. Верхний узел, из которого начинается частотное разделение, имеет нулевой номер. Последний уровень имеет номер Nw, означающий число делений частотных диапазонов. Таким образом, на каждом уровне, в том числе и на последнем, ряд узлов можно представить как последовательность вейвлет-спектров, спектры которых разделяют весь частотный диапазон от нуля до частоты Найквиста.
Узлы обозначаются следующим образом [Nw, n], где Nw – уровень декомпозиции, n – порядковый номер узла на заданном уровне Nw.
Ожидается, что спектры вейвлет-коэффициентов узлов идут по частоте один за другим по мере увеличения порядкового номера узла.
В данном случае прореженные вейвлет-коэффициенты подвергаются процедуре реконструкции для восстановления числа отсчетов. Это делается для того, чтобы спектры вейвлет-коэффициентов соответствовали частотному диапазону исходного сигнала. Тогда спектр вейвлет-коэффициентов сигнала начального узла [Nw, 0] будет находиться в начале частотного диапазона, а спектр вейвлет-коэффициентов сигнала конечного узла [Nw, 2^Nw-1] будет находиться в конце частотного диапазона, примыкать к частоте 5 кГц (в данном примере частота дискретизации Fs равна 10 кГц).
Ниже приводится код, который позволяет получить диаграмму спектров в осях: ордината – номера узлов, абсцисса – частота спектров. Уровень декомпозиции Nw=5.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
echo on
clc
echo off
%% Сигнал
Fs=10000; % Частота дискретизации
S=wgn(1, 50*Fs, 0); % Белый шум
N=length(S); % Длина вектора сигнала
Ffft=0:Fs/N:Fs-Fs/N; % Частота для спектра Фурье
%% WP
Nw=5; %% Уровень декомпозиции
wavelet = 'db40'; % Вейвлет
T = wptree(2,Nw,S,wavelet); % Дерево декомпозиции пакетного вейвлета
%% Цикл для вычисления схемы частотного положения спектров
% вейвлет-коэффициентов на последнем уровне Nw=5
for i=0:1:2^Nw-1 %% i - порядковый номер узла на заданном уровне Nw
rcfs = wprcoef(T,[Nw i]); % Вейвлет-коэффициенты в заданном узле i
Yrcfs=abs(fft(rcfs)); % Спектр фурье вейвлет-коэффициентов
% Так как спектры оказались сгруппированными по четыре, то решено
% выполнить их окраску «по четверкам»
% Цвета спектров каждой четверки:
% ЧЕРНЫЙ – узел с номером [5, 0], т.е. с номером 0,
% а также каждый следующий узел с номером «плюс 4»
% ГОЛУБОЙ – узел с номером 1 и далее – каждый с номером «плюс 4»
% ЗЕЛЕНЫЙ – узел с номером 2 и далее – каждый с номером «плюс 4»
% КРАСНЫЙ – узел с номером 3 и далее – каждый с номером «плюс 4»
% Цветовая дифференциация позволяет обратить внимание на инверсию
% порядка следования третьего (зеленый) и четвертого (красный)
% спектра внутри каждой четверки.
% И не только, имеет место инверсия всей четверки.
% Зачем, не понятно
vect32=['k' 'b' 'g' 'r' ...
'k' 'b' 'g' 'r' ...
'k' 'b' 'g' 'r' ...
'k' 'b' 'g' 'r' ...
'k' 'b' 'g' 'r' ...
'k' 'b' 'g' 'r' ...
'k' 'b' 'g' 'r' ...
'k' 'b' 'g' 'r'];
vect=[vect32 vect32 vect32 vect32];
%% PLOT
plot(Ffft, i+Yrcfs./max(Yrcfs), vect(i+1)); grid on; hold on; title('Nw=5; wavelet=db40');
xlabel('F, Hz'); ylabel('Node number');
set(gca,'FontName','TimesET'); set(gca,'FontWeight', 'normal');set(gca,'FontSize',14);
axis([0 5000 0 i+1]);
end
%% END %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Рисунок, полученный с помощью этого кода, размещен по адресу:
http://www.imageup.ru/img2/3408073/wp_w ... y.jpg.html
На рисунке над спектрами для улучшения чтения дополнительно проставлены номера узлов, которые соответствуют номерам на оси ординат.
Из рисунка видно, что спектры узлов не располагаются по оси частот последовательно, один за другим, сообразно номерам их узлов. Самым крайним спектром (справа, вблизи частоты 5 кГц) оказался спектр узла с номером 16, а спектр последнего узла (с номером 31) оказался в средине частотного диапазона – 3300-3500 Гц.
Если исходить из трактовки усовершенствованного пирамидального алгоритма Мала (Суть алгоритма состоит в разбиении частотного диапазона сигнала на две части: высокочастотную и низкочастотную, каждая из которых в свою очередь также разбивается на высокочастотную и низкочастотную части и так – до достижения заданного уровня), то частотные диапазоны спектров узлов должны идти один за другим по мере увеличения порядкового номера узла.
Однако, все перепутано местами.
Коллеги, подскажите, пожалуйста, как добиться, чтобы спектры вейвлет-коэффициентов узлов шли по частоте один за другим соответственно номерам их узлов.
С уважением, Вениамин.
Re: Вейвлет-спектр энергии шума имеет ступени
как добиться, чтобы спектры вейвлет-коэффициентов узлов шли по частоте один за другим соответственно номерам их узлов.
Просто переставьте их в нужном порядке. Последовательность номеров на вашем рисунке образует код Грея. Получить таблицу преобразования можно, например, так:
Код: Выделить всё
ind = bin2gray(0:31, 'pam', 32)
ind =
Columns 1 through 13
0 1 3 2 6 7 5 4 12 13 15 14 10
Columns 14 through 26
11 9 8 24 25 27 26 30 31 29 28 20 21
Columns 27 through 32
23 22 18 19 17 16
С уважением
Александр Сергиенко
Александр Сергиенко
Re: Вейвлет-спектр энергии шума имеет ступени
Спасибо, Александр Борисович, за консультацию.
Все же хотелось бы понять, как из стройного алгоритма деления частотных диапазонов низкочастотным L_o и высокочастотным H_o фильтрами получилась путаница узлов.
С уважением, Вениамин
Все же хотелось бы понять, как из стройного алгоритма деления частотных диапазонов низкочастотным L_o и высокочастотным H_o фильтрами получилась путаница узлов.
С уважением, Вениамин
Re: Вейвлет-спектр энергии шума имеет ступени
Не забывайте, что после деления частотного диапазона на две части происходит и понижение частоты дискретизации. Представьте, что при этом произойдет с частотами в высокочастотной половинке (с учетом неоднозначности понятия частоты для дискретного сигнала).
С уважением
Александр Сергиенко
Александр Сергиенко
Re: Вейвлет-спектр энергии шума имеет ступени
Да, я помню. Для адекватного восприятия спектров вейвлет-коэффициентов выполнена реконструкция децимированных вейвлет-коэффициентов. В этом случае спектр исходного сигнала и спектры вейвлет-коэффициентов можно адекватно сравнивать на одной оси частот от 0 до частоты Найквиста.
Возможно, где-то в коде дерева вейвлет-преобразования включена процедура преобразования в код Грея. По-видимому для улучшения процедуры сжатия информации.
Спасибо, Александр Борисович, за консультацию.
С уважением, Вениамин.
Возможно, где-то в коде дерева вейвлет-преобразования включена процедура преобразования в код Грея. По-видимому для улучшения процедуры сжатия информации.
Спасибо, Александр Борисович, за консультацию.
С уважением, Вениамин.
