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

[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

[fontp]

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

Слишком много слов. Чтобы иметь РАЗРЕШЕНИЕ в классическом смысле Вы должны взять 2 гармоники равной амплитуды неразрешимые исходным ДПФ (действительно разность частот 0.5 в Ваших обозначениях), причём некратных шагу ДПФ. Берите +0.33 и +0.83 для начала. Потом разрешать Вы их должны всегда (с вероятностью 1, любые частоты f1 и f1+0.5), поскольку автору нужно иметь измерительный прибор, а не пробник.

Вот эту ерунду Вы обещали получить. Но нарисовали единственную гармонику.

Остальное, что я написал про фазовые шумы - это попытка дать Вам шанс, упростив задачу. Забыть о РАЗРЕШЕНИИ (оно ведь действительно не нужно автору так же как и точность измерения часоты), а попытаться измерять расширение спектральной линии ВЫЗВАННОЕ фазовыми шумами у единственной гармоники. При определённой модели фазовых шумов, может быть, что-то и получится. Но как инструментальное измерительное средство - никак, пробничек получится.

Пример: аналоговый генератор у которого фаза прыгает на ПИ, как обычно, но не в случайные моменты времени, а периодически, примерно так же как и у измерителя. Фишка так легла, что все скачки фазы легли между блоками измерений. Ваш прибор покажет, что генератор идеален и фазовые шумы отсутствуют, поскольку расширение спектральной линии соответствуют погрешности прибора.

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

Пост отредактирован. fontp - устное предупреждение за хамство. Еще одно такое полуоскорбление - будет readonly. DS

[Stanislav]

Слишком много слов.

Это у меня-то?
Я пишу здесь по делу, и слов употребляю ровно столько, сколько нужно для доходчивого пояснения.
Не моя вина, что некоторым объяснять нужно гораздо более одного раза.
Лучше обратите внимание на сообщения некоторых других участников темы. Вот где уж действительно "вольница мысли"...
И ещё: не поленитесь, пожалуйста, уделить внимание изучению английского языка. В части произношения имён собственных шотландского происхождения. Возможно, желания писать глупости после этого поубавится.

...Чтобы иметь РАЗРЕШЕНИЕ в классическом смысле Вы должны взять 2 гармоники равной амплитуды неразрешимые исходным ДПФ (действительно разность частот 0.5 в Ваших обозначениях), причём некратных шагу ДПФ.

Пожалуйста, не нужно адресовать мне, да ещё кричащими буквами, прописную истину, которую я безуспешно пытался пояснить Вам же с самого начала темы (см. посты №№ 6,8 и т.д, вплоть до цитируемого Вами). Если считаете, что понимание всё-таки наступило, пожалуйста, не нужно выражать его столь эмоционально.
Я же так не считаю. Посему, "ликбез" придётся продолжить... Дальше прошу читать внимательно: повторяться более не желаю.

Если какая-либо линейная система может отделить гармонический сигнал определённой амплитуды от другого гармонического сигнала гораздо бОльшей амплитуды, это означает, что её разрешающая способность по частоте не является худшей разности истинных частот данных сигналов. ©Stanislav.

Возражения имеются? Если да - только на примере: словесный спам достал уже.

...Берите +0.33 и +0.83 для начала. Потом разрешать Вы их должны всегда (с вероятностью 1, любые частоты f1 и f1+0.5), поскольку автору нужно иметь измерительный прибор, а не пробник.

Вот эту ерунду Вы обещали получить. Но нарисовали единственную гармонику.

Во-первых, я никому, в том числе, и Вам, разрешать ничего не должен.
Во-вторых, как это "единственную", если их там две?
Ей-богу, как дитё малое... Сначала дай две частоты, а теперь ещё и третью... Эдак не напасёшься...
Могли бы и сами это проделать. Я не прошу изучать внимательно программу, хотя там и каменты в каждой строчке, но "сгенерить" ещё одну синусоиду через копи-паст, наверное, можно было и самостоятельно.

Так уж и быть, проделаю это за Вас, c названным отклонением частот. Интервал выбран <1мс.


- стык кусков до коррекции фазы основной частоты,

- стык кусков после коррекции фазы основной частоты,

- модуль спектральной функции одного куска,


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

...Остальное, что я написал про фазовые шумы - это попытка дать Вам шанс, упростив задачу.

Да бросьте.
Не нужно мне "давать шанс", равно как и "упрощать задачу". Тем более, что это не моя задача, а задача автора темы. Так что прошу обратить внимание на задачу, и впредь писать по существу... кхм... дела.

...Забыть о РАЗРЕШЕНИИ (оно ведь действительно не нужно автору так же как и точность измерения часоты),

Да что Вы говорите? Маклауда, стало быть, уже побоку?

Опять сильное заявление. К которому нужно, как я уже и объяснял, прибавлять "по моему мнению".
Читайте внимательнее постановку задачи. Можно начать с названия темы.

ЗЫ. Всеобъемлемость Вашего утверждения легко опровергнуть примером. Но - по запросу.

...а попытаться измерять расширение спектральной линии ВЫЗВАННОЕ фазовыми шумами у единственной гармоники. При определённой модели фазовых шумов, может быть, что-то и получится. Но как инструментальное измерительное средство - никак, пробничек получится.

Пример: аналоговый генератор у которого фаза прыгает на ПИ, как обычно, но не в случайные моменты времени, а периодически, примерно так же как и у измерителя. Фишка так легла, что все скачки фазы легли между блоками измерений. Ваш прибор покажет, что генератор идеален и фазовые шумы отсутствуют, поскольку расширение спектральной линии соответствуют погрешности прибора.

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

Перестаньте кривляться

Не пытался, даже ничуть. Отвечал на конкретный тезис известного барона по поводу "экстраполяции фазы".

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

Прошу прощенья за оффтоп...

[GetSmart]

- стык кусков до коррекции фазы основной частоты,

- стык кусков после коррекции фазы основной частоты,

- модуль спектральной функции одного куска,


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

И что мы видим? В исходных данных было две разных, но близких частоты (примем на веру) одинаковой амплитуды. На спектре склейки только одна явная вершина. Остальные два пика маленькие, грубо ~50 раз меньше основного. Почему их 2 и почему в 50 раз меньше - хз. Вобщем окончательный спектр представляет совсем не то, что нужно автору. Кроме того, этот спектр до неприличия похож на спектр склейки одинаковых синусоид. Да и частоты этих пиков на спектрограмме слишком отдалены друг от друга.