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

[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. Или шум пропущенный через узкополосный фильтр. Нужно, оценить уширение спектра с целью определить параметры стабильности. Что-то в этом духе

[GetSmart]

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

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

А главное, во всём виноват Гильберт с его заумным фильтром

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

Какие параметры (качества) сигнала Вы хотите получить от метода измерения? И до сих пор не ясно есть ли в сигнале кроме основной синусоиды другие сигналы? (ниже -60 дб не важно)

По поводу "как склеить М кусков чтобы получить дополнительную информацию". Вас не устроит в качестве результата статистика максимально точной частоты основной грамоники в каждом из кусков и её среднеквадратичного отклонения? Чем больше будет накоплено кусков, тем более точный будет результат вычисления средней частоты (основного пика в FFT). Если не пропорционально кол-ву кусков, то хотя бы корню квадратному из их кол-ва. Причём, здесь уже не будет проблемы с ограниченностью псевдокоггерентного FFT, у которого точность никогда не превысит точность склеек, которая может быть намного хуже разрешающей способности FFT M*N.

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