Увеличение разрешения по частоте, засчёт накопления Опции

[EKirshin]

Здравствуйте!

Такая задачка.

Есть стационарный сигнал: шум + синусоида.

Из этого сигнала периодически набирается N отсчётов.
Вопрос: можно ли, имея M кусков этого сигнала использовать их для получения спектра сигнала размерностью M*N?

Какой будет эффект, если на каждый из этих кусков наложить окно, "склеить" эти куски (получится сигнал длиной M*N), затем сделать fft на M*N точек?

[Stanislav]

Вот обещанная "рыба".

Код

function phase_adj()
fsampl = 1000;      % частота дискретизации, кГц
freq = 68;          % Частота гармонического тона, кГц
len_buf = 4096;     % длина буфера, отсчёты
shift = pi/4;       % величина фазового сдвига между кусками, рад

% Заполняем буферы сигналом и аддитивным шумом
buffer_1 = sin(2*pi*(1:len_buf)*freq/fsampl)+(1e-5)*randn(1,len_buf);
buffer_2 = sin(2*pi*(1:len_buf)*freq/fsampl+shift)+(1e-5)*randn(1,len_buf);

% Пытаемся склеить без учёта фазовых соотношений
sum_buf_1 = [buffer_1, buffer_2];

% Делаем оценку спектров кусков и склейки, используя спектральные окна
specw_1=fft(blackman(len_buf)'.*buffer_1);
specw_2=fft(blackman(len_buf)'.*buffer_2);

sp_sum_1=fft(hamming(2*len_buf)'.*sum_buf_1);

% Находим фазовые спектры кусков
ang_1 = angle(specw_1);
ang_2 = angle(specw_2);

% Находим положение максимума модуля спектральной функции
[max_1, frq_1] = max(abs(specw_1));

% Находим фазу основного тона в конце 1-го куска (не доделано!)
ang_1_end = ang_1(frq_1);

% Находим разность фаз
diff_ang=ang_1_end-ang_2(frq_1);

% Находим поворачивающий множитель для частоты основного тона
rotator=-exp(i*diff_ang);

% Преобразуем второй буфер в комплексный вид
buf_2_hilb=hilbert(buffer_2);

% Разворачиваем его на величину сдвига фаз (нужно ещё ввести компенсацию по
% частоте, проще говоря - временнОй сдвиг)
buf_2_rot=real(buf_2_hilb*rotator);

% Склеиваем
sum_buf_2 = [buffer_1, buf_2_rot];

% Находим модуль спектра склейки
sp_sum_2=fft(hamming(2*len_buf)'.*sum_buf_2);

% Выводим результаты
figure (1)
plot (buffer_1(1:100))
hold on
plot (buffer_2(1:100), 'r')
hold off
grid on
title ('Signals')

figure (2)
plot ((4032:4160), sum_buf_1(4032:4160))
grid on
title ('Glue 1')

figure (3)
plot ((4032:4160), sum_buf_2(4032:4160), 'Color', [0.9 0 0])
grid on
title ('Glue 2')

figure (4)
semilogy (500:600, abs(sp_sum_1(500:600)))
hold on
semilogy (500:600, abs(sp_sum_2(500:600)), 'r')
hold off
grid on
title ('Spectra')


return;

Правда, она годится лишь для демонстрации принципа: для условий задачи она пока не слишком хорошо подходит - не сделан временнОй сдвиг, поэтому, компоненты, отстоящие от основной частоты на значительную величину, будут "склеиваться" плохо. Кроме того, фаза определяется весьма грубо, и нет процедуры её коррекции от перескока. Однако, вблизи центральной частоты всё будет более-менее хорошо.
За ошибки прошу не карать строго - набросал относительно быстро, и они возможны.

Завтра постараюсь сделать аккуратно, со всеми необходимыми компенсациями. Скорее всего - в спектральной области.
Ниже приведены картинки.

- фрагменты исходных сигналов.

- склейка без разворота фазы. Виден её скачок.

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


- фрагмент спектров склеек. Синий - без разворота, красный - с разворотом. Разрешение по частоте увеличилось по сравнению с исходными кусками в 2 раза!

Сообщение отредактировал Stanislav - Aug 2 2008, 07:52

[fontp]

- фрагмент спектров склеек. Синий - без разворота, красный - с разворотом. Разрешение по частоте увеличилось по сравнению с исходными кусками в 2 раза!

Маклауд такие синусоиды считает значительно точнее, в разы.

Вы же сами утвеждали, что разрешение для идеальной одиночной синусоиды не определено.
Дайте две.

Или лучше внесём в модель фазовые шумы, возьмём толстый фломастер и напишем на крышке прибора:

" Прибор для измерения фазовых шумов генератора посредством экстраполяции фазы"
и подпись

"барон Мюнхаузен"

ЗЫ. Вообще-то постановка задачи у автора опять не получилась. Если идеальная синусоида - то спектр её делта-функция. Осталось оценить частоту. Как было сказано - это не тот случай.
Значит скорее всего - это узкополосный сигнал, например сигнал аналогового генератора. Или DDS. Или шум пропущенный через узкополосный фильтр. Нужно, оценить уширение спектра с целью определить параметры стабильности. Что-то в этом духе

[Stanislav]

Маклауд такие синусоиды считает значительно точнее, в разы.

Странно, что Вы упорно не желаете понять простой вещи: Автору темы нужно не точное измерение частоты гармонического сигнала, а получение модуля спектральной функции с разрешением по частоте более высоким, чем даёт каждая из накопленных последовательностей.
Найти частоту гармонического сигнала с хорошей точностью при таком отношении С/Ш можно и без Маклауда.
ЗЫ. Так Маклауд, или всё-таки "не-Маклауд"? Как по-русски писать правильно?

...Вы же сами утвеждали, что разрешение для идеальной одиночной синусоиды не определено.
Дайте две.

Это был только набросок алгоритма, ещё не работающий, как положено. Две частоты дам ниже.

...Или лучше внесём в модель фазовые шумы, возьмём толстый фломастер и напишем на крышке прибора:

" Прибор для измерения фазовых шумов генератора посредством экстраполяции фазы"
и подпись

"барон Мюнхаузен"

Чушь.
Шумы в модели присутствуют. Для проверки разрешения синусоид достаточно.
Кроме того, Вы не поняли сути предложенного метода: ни о какой "экстраполяции фазы" речь здесь не идёт.

...ЗЫ. Вообще-то постановка задачи у автора опять не получилась. Если идеальная синусоида - то спектр её делта-функция. Осталось оценить частоту. Как было сказано - это не тот случай.

Что Вам не нравится в постановке? По-моему, теперь все необходимые данные есть.

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

Ну, наконец-то...
Я об этом с начала темы талдычу. Да и EKrishin писал недвусмысленно, что ему нужен именно спектр сигнала вблизи основной частоты. Просто удивительно, что там могло быть непонятного?

.............................................................................
Вот доработанный алгоритм. С ним уже можно поиграться - попробовать менять частоту и сдвиг фаз. Частота предполагается известной, чтобы не загромождать код. Второй гармонический тон меньшей амплитуды добавлен.

Код

function phase_adj()  
fsampl = 1000;      % частота дискретизации, кГц
freq = 40;          % частота гармонического тона, кГц

len_buf = 4096;     % длина буфера, отсчёты
shift = pi;       % величина фазового сдвига между кусками, рад

% Заполняем буферы сигналом и аддитивным шумом
buffer_1 = sin(2*pi*(0:len_buf-1)*freq/fsampl)+(1e-5)*randn(1,len_buf);
buffer_1 = buffer_1+0.03*sin(2*pi*(0:len_buf-1)*(freq+0.5)/fsampl);

shift_1=shift; % не корректно; токмо для ускорения процесса

buffer_2 = sin(2*pi*(0:len_buf-1)*freq/fsampl+shift)+(1e-5)*randn(1,len_buf);
buffer_2 = buffer_2+0.03*sin(2*pi*(0:len_buf-1)*(freq+0.5)/fsampl+shift_1);

% Пытаемся склеить без учёта фазовых соотношений
sum_buf_1 = [buffer_1, buffer_2];

% Делаем оценку спектров кусков и склейки, используя спектральные окна
specw_1=fft(hamming(len_buf)'.*buffer_1);
specw_2=fft(hamming(len_buf)'.*buffer_2);
sp_sum_1=fft(hamming(2*len_buf)'.*sum_buf_1);

% Находим фазовые спектры кусков
ang_1 = angle(specw_1);
ang_2 = angle(specw_2);

% Находим положение максимума модуля спектральной функции
[max_1, frq_1] = max(abs(specw_1));

% Находим разность фаз
diff=ang_1(frq_1)-ang_2(frq_1);
dif_ang=mod(diff+2*pi*freq/fsampl*len_buf, 2*pi);

% Находим поворачивающий множитель для частоты основного тона
rotator=exp(i*dif_ang);

% Преобразуем второй буфер в комплексный вид
buf_2_hilb=hilbert(buffer_2);

% Разворачиваем его на величину сдвига фаз (более правильно делать сдвиг во времени)
buf_2_rot=real(buf_2_hilb*rotator);

% Склеиваем
sum_buf_2 = [buffer_1, buf_2_rot];

% Находим модуль спектра склейки
sp_sum_2=fft(hamming(2*len_buf)'.*sum_buf_2);

% Выводим результаты
figure (1)
plot (1:100, buffer_1(1:100))
hold on
plot (1:100, buffer_2(1:100), 'r')
hold off
ylim([-1.5 1.5]);
grid on
title ('Signals')

figure (2)
plot ((4046:4146), sum_buf_1(4046:4146))
ylim([-1.5 1.5]);
grid on
title ('Glue 1')

figure (3)
plot ((4046:4146), sum_buf_2(4046:4146), 'Color', [0.9 0 0])
ylim([-1.5 1.5]);
grid on
title ('Glue 2')

figure (4)
semilogy (1:1000, abs(sp_sum_2(1:1000)))
grid on
title ('Spectrum of long sequence')

figure (5)
semilogy (1:1000, abs(specw_1(1:1000)))
grid on
title ('Spectrum of short sequence')
return;

fontp, наблюдайте внимательно:

- модуль спектра последовательности одинарной длины,


- модуль спектра последовательности двойной длины, полученной путём "склеивания".

И никаких Маклаудов.


...................................................
Далее, если позволит время, попробую сделать то же самое через Фурье - должно получиться красивее.

Сообщение отредактировал Stanislav - Aug 2 2008, 17:03

[alex_os]

- модуль спектра последовательности двойной длины, полученной путём "склеивания".

И никаких Маклаудов.
...................................................
Далее, если позволит время, попробую сделать то же самое через Фурье - должно получиться красивее.

Фаза второй палки (shift_1) у Вас удачно попала . Если взять shift_1 = 1 то от второй палки рожки да ножки остаются ...

[Stanislav]

Фаза второй палки (shift_1) у Вас удачно попала . Если взять shift_1 = 1 то от второй палки рожки да ножки остаются ...

Это верно. Там сдвиг фаз второй синусоиды берётся "от балды", о чём и написано в каменте (см. текст программы), и поэтому подобрана "хорошая" частота второй синусоиды. Но, если посчитать её фазу аккуратно, а потом ещё ввести при формировании комплексного сигнала коррекцию фазы по частоте (или, проще, сдвиг во времени), то всё будет в ажуре. Корректирующий фильтр считать лень, но, если время будет, сделаю... Попробую в базисе Фурье алгоритм реализовать, только попожже.
На коротких кусках же разрешить вторую частоту вообще не удаётся...

[alex_os]

Это верно. Там сдвиг фаз второй синусоиды берётся "от балды", о чём и написано в каменте (см. текст программы), и поэтому подобрана "хорошая" частота второй синусоиды. Но, если посчитать её фазу аккуратно, а потом ещё ввести при формировании комплексного сигнала коррекцию фазы по частоте (или, что проще, сдвиг во времени), то всё будет в ажуре. Корректирующий фильтр считать лень, но, если время будет, сделаю... Попробую в базисе Фурье алгоритм реализовать, только попожже.
На коротких кусках же разрешить вторую частоту вообще не удаётся...

А откуда Вы узнаете сдвиг во времени если задержка между блоками не известна?!

[Stanislav]

А откуда Вы узнаете сдвиг во времени если задержка между блоками не известна?!

Измерить!
Ну, ей-богу... Если какие-то величины не известны, кто мешает их получить из сигнала, или самого процесса ввода данных? Такая простая мысль, как видно, не посещает многих читателей умных книжек.
Как я и писал уже, слишком точно временнОй сдвиг измерять и не нужно. Но детали следует обдумать.