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

[EKirshin]

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

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

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

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

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

[alexkok]

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

Я думаю что это не реализуемо.
Метод, который предлагает Stanislav - подкручивание фазы основной составляющей, будет работать (возможно, не уверен) только для составляющих коррелированных с основной гармоникой.
Для всех прочих он работать не будет, т.к. они будут складываться с произвольными фазами.
Если есть возможность получать время начала каждой пачки с точностью до отсчёта, то можно вставить пустые (забитые нулями) промежутки и расфильтровать суммарную выборку.
На каждый из этих "N" кусков наложить окно перед склейкой.

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

[Stanislav]

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

Это основное ограничение, никуда не деться. Однако, остальное можно и нужно считать шумом, а как мы выяснили, он эргодичен.

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

Код

function phase_adj()
fsampl = 1000;      % частота дискретизации, кГц
freq = 41;          % Частота гармонического тона, кГц
len_buf = 4096;     % длина буфера, отсчёты
gap = 0.35;          % интервал между кусками, мс, не более 8
ss_fact = 5;        % количество точек модели сигнала на 1 выборку АЦП

% Формируем опорную сетку для "реального" сигнала
tick = ss_fact*len_buf;

% Формируем модель "реального" сигнала с разрешением бОльшим, чем выборка АЦП
% (коряво: можно получить куски и аналитически, но это для сравнения спектров)
signal = sin(2*pi*(0:4*tick)*freq/(ss_fact*fsampl)) + (1e-2)*randn(1, 4*tick+1) +...
    0.02*sin(2*pi*(0:4*tick)*(freq+0.5)/(ss_fact*fsampl)) +...
    0.02*sin(2*pi*(0:4*tick)*(freq-0.5)/(ss_fact*fsampl));

% Заполняем буферы
buffer_1 = signal(1:ss_fact:tick);
sampl_gap = floor(1e6*ss_fact*gap/fsampl);
buffer_2 = signal(tick+1+sampl_gap:ss_fact:tick+sampl_gap+tick);

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

% Делаем оценку спектров кусков и "реального" сигнала, используя спектральные окна
specw_1 = fft(chebwin(len_buf, 60)'.*buffer_1);
specw_2 = fft(chebwin(len_buf, 60)'.*buffer_2);
sp_sig =  fft(chebwin(2*len_buf, 60)'.*signal(1:ss_fact:2*tick));

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

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

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

% Находим временнУю задержку, соответствующую разности фаз
tau = dif_ang/(2*pi*freq);

% Находим величину сдвига в отсчётах сигнала
shift = round(tau*fsampl);

% Сдвигаем 2-й кусок для слияния фаз
tmp_buff = buffer_2(shift+1:end); % весьма грубо, но пока что сойдёт

% Kоличество отсчётов на период (грубо)
period = round(fsampl/freq);

% Дописываем "хвост"
tmp_buff = [tmp_buff, tmp_buff(end-period+1:end-period+shift)]; % пуркуа бы и не па?

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

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

% Выводим результаты
bins = (0:len_buf-1)*fsampl/len_buf;

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 (bins(1:1000), abs(specw_1(1:1000)))
ylim([0.1 1000]);
grid on
xlabel('Frequency, Hz')
title('Spectrum of short sequence')

figure (5)
semilogy (bins(1:1000)/2, abs(sp_sig(1:1000)), 'b')
hold on
semilogy (bins(1:1000)/2, abs(sp_sum_2(1:1000)), 'r')
hold off
ylim([0.1 1000]);
grid on
xlabel('Frequency, Hz')
legend('Real signal', 'Adjusted pieces')
title('Spectra of long sequences')
return;

Пояснение: сначала формируется сигнал с временнЫм разрешением 1/5 относительно частоты выборки (имитация "непрерывнрсти"), а затем из него формируются исходные последовательности, с учётом задержки между блоками (gap). Стыковка производится с точностью в +/- 0,5 длины отсчёта. Повысить её можно путём интерполяции.

...Для всех прочих он работать не будет, т.к. они будут складываться с произвольными фазами.
Если есть возможность получать время начала каждой пачки с точностью до отсчёта, то можно вставить пустые (забитые нулями) промежутки и расфильтровать суммарную выборку...

Мне кажется, что измерять временнОй интервал с точностью до фазы основной гармоники вовсе и не нужно - требуется лишь получить погрешность не более ~ 1/(2*F ), где F - ширина интересующего диапазона.
Небольшие нестыковки фаз компонентов сигнала на результат качественно не повлияют; будет только небольшая неопределённость амплитуд и частот спектральных составляющих. Но разрешение как таковое будет увеличено кардинально.

...Если же время начала получить невозможно, то можно лишь слегка улучшить разрешение если до БПФ вычесть из сигнала основную составляющую и фильтровать через БПФ только остаток.
Такой метод применяется при аналоговых измерениях чистоты спектра. Преимуществом при этом будет возможность применения более узкополосного окна или БПФ без окна вообще...

Да, это верно, но всё-таки не принципиально. Для подобных измерений подавление тестового сигнала обячно делается из-за недостаточной разрешающей способности аппаратуры по уровню (или, более правильно, недостаточного динамического диапазона).
Имея в распоряжении, скажем, АЦП со 120-130 дБ ДД, подавлять основную частоту, скорее всего, и не потребуется.
А при узкополосной фильтрации кусок сигнала просто "развалится" по причине большой добротности фильтра, и никакого выигрыша, скорее всего, не получится...

...Если есть априорная информация об измеряемых составлящих - как правило это гармоники сети и частоты сравнения в синтезаторах, то можно ещё попытаться улучшить измерение, но не БПФом.

Совершенно верно, я об этом тоже думал.
Если сигнал имеет какие-либо устойчивые артефакты помимо основного тона, делу разрешения неодозначностей это может сильно помочь, и даже без измерения интервала времени вообще - экстремальный подход к проблеме. Однако, они должны лежать заметно выше шумового фона.
Более того, уместно было бы ещё хотя б один пилотный тон известной частоты во входной сигнал добавить. Я бы, пожалуй, выбрал эту частоту в районе 100-500Гц.
Однако, и простейший счётчик интервалов времени, думается, работать будет очень хорошо. Попробую ввести в модель пожже.

Собственно, ПФ не является чем-то необходимым. Требуемые преобразования можно осуществить целиком и во временнОй области (методами регрессии, например), которые ничем не хуже ДПФ.

Вот картинки. Хорошо видно, что боковые лепестки (+/- 0,5 кГц) в склейке разрешаются так же хорошо, как и в исходном сигнале, при условии, что последний анализируется на двойной длине куска. В то время, как на куске одинарной длины "боковушки" не разрешаются вообще.
Алгоритм работает практически с любыми настройками, но, ввиду грубости стыковки кусков, частоту "высокой" делать не стОит.