Перейти к содержанию
    

Увеличение разрешения по частоте

Дык если можете делать фильтрацию - децимацию, так уж фильтруйте и децимируйте до полосы скажем 20 кГц и соотв. частоте дискретизации. И на низкой Fs берете 4К отсчетов и анализируете спокойненько. (Зачем вам Fs = 1 МГц если полезная полоса 10КГц ?!).
А какой смысл? Разрешения по частоте это не увеличит; можно будет только вычислительный ресурс сэкономить, но в контексте данной темы это не принципиально.

 

...Попытка сшивания и фазирования блоков приведет скорее всего к паразитной фазовой модуляции в сшитом блоке, т.е. интересующая Вас палка в спектре размажется...
А вот и посмотрим, как она размажется. :)

Смысл "сшивки" как раз и состоит в том, чтобы устранить паразитные перескоки фаз между блоками.

 

...Даже более того преставим что рядом с основной гармоникой есть "палка" меньшего уровня , когда сошьем блоки так чтобы основная гармоника была без разрыва фазы, "палка" же будет иметь разрывы фазы, в результате в сшивке вместо истинного спектра увидим основную гармонику и спектр фазоманипулированной палки , те ХЗЧ.
А если хорошо подумать? ;)

Близкие к основной "палке" компоненты спектра когерентны по отношению к ней - сигнал стационарен. Пересчёт поворачивающих множителей для них не составит труда. :)

 

...Или пойти еще дальше возмем одну выборку (будем считать что она без шума) и n раз сшив ее саму с собой получим здоровую сшивку :) - увидим ли в спектре этой хрени чего-то новое? Максимум что получится это более точно измерить частоту основной гармоники но это можно сделать и гораздо более простыми методами.
Ну, это уже лишнее. :biggrin:

...Так что природу не обманешь.
Это точно.

Поделиться сообщением


Ссылка на сообщение
Поделиться на другие сайты

А какой смысл? Разрешения по частоте это не увеличит; можно будет только вычислительный ресурс сэкономить, но в контексте данной темы это не принципиально.

Дык я думал что разрешение по частоте пропорционально Fs/Nfft :).

Поделиться сообщением


Ссылка на сообщение
Поделиться на другие сайты

Дык я думал что разрешение по частоте пропорционально Fs/Nfft :) .
Дык, Вы ж децимировать предлагаете. :) Сколько выборок останется от 4096?

 

ЗЫ. Понял. Вы предлагаете брать все 4096 отсчётов на НЧ, увеличив время наблюдения. Что ж, тоже вариант. Не знаю, однако, насколько реализуемый и удобный для Автора темы.

 

Лана, пора делом заняться... :biggrin:

 

ЗЗЫ. Посмотрел файлы. Похоже, что у Автора темы всё именно так и сделано - сигнал очень низкочастотный, и в квадратурах.

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

2 EKrishin

Скажите, в файлах, выложенных Вами, частота выходных отсчётов действительно 1 МГц? Или всё-таки менее?

Изменено пользователем 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;

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

За ошибки прошу не карать строго - набросал относительно быстро, и они возможны. :)

 

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

Ниже приведены картинки.

 

post-4987-1217634850_thumb.jpg - фрагменты исходных сигналов.

 

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

 

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

 

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

Изменено пользователем Stanislav

Поделиться сообщением


Ссылка на сообщение
Поделиться на другие сайты

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

 

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

 

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

Дайте две.

 

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

 

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

и подпись

 

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

 

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

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

Поделиться сообщением


Ссылка на сообщение
Поделиться на другие сайты

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

и подпись

 

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

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

 

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

...

Под стационарностью я понимаю неизменность параметров сигнала: частота сигнала (синусоиды) остаётся неизменной (не учитывая качество самого сигнала: дрожание частоты и т.п., изменяющие частоту сигнала), уровень шума и его другие характеристики также неизменны.

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

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

 

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

Изменено пользователем GetSmart

Поделиться сообщением


Ссылка на сообщение
Поделиться на другие сайты

Маклауд такие синусоиды считает значительно точнее, в разы.
Странно, что Вы упорно не желаете понять простой вещи: Автору темы нужно не точное измерение частоты гармонического сигнала, а получение модуля спектральной функции с разрешением по частоте более высоким, чем даёт каждая из накопленных последовательностей.

Найти частоту гармонического сигнала с хорошей точностью при таком отношении С/Ш можно и без Маклауда.

ЗЫ. Так Маклауд, или всё-таки "не-Маклауд"? Как по-русски писать правильно? ;)

 

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

Дайте две.

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

 

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

 

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

и подпись

 

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

Чушь.

Шумы в модели присутствуют. Для проверки разрешения синусоид достаточно.

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

 

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

...Значит скорее всего - это узкополосный сигнал, например сигнал аналогового генератора. Или 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, наблюдайте внимательно:

 

post-4987-1217689806_thumb.jpg- модуль спектра последовательности одинарной длины,

 

post-4987-1217689818_thumb.jpg- модуль спектра последовательности двойной длины, полученной путём "склеивания". ;)

 

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

 

 

...................................................

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

Изменено пользователем Stanislav

Поделиться сообщением


Ссылка на сообщение
Поделиться на другие сайты

post-4987-1217689818_thumb.jpg- модуль спектра последовательности двойной длины, полученной путём "склеивания". ;)

 

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

...................................................

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

 

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

 

post-17030-1217697363_thumb.jpg

Поделиться сообщением


Ссылка на сообщение
Поделиться на другие сайты

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

На коротких кусках же разрешить вторую частоту вообще не удаётся...

Поделиться сообщением


Ссылка на сообщение
Поделиться на другие сайты

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

На коротких кусках же разрешить вторую частоту вообще не удаётся...

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

Поделиться сообщением


Ссылка на сообщение
Поделиться на другие сайты

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

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

Как я и писал уже, слишком точно временнОй сдвиг измерять и не нужно. Но детали следует обдумать.

Поделиться сообщением


Ссылка на сообщение
Поделиться на другие сайты

Используем более быстрые АЦП, затем производится фильтрация с децимацией, в результате чего получаем повышенное SNR.

 

Длительность записи - 4096 точек на 1 Мгц. Длительность паузы неопределена. Порядок - около 100 мс.

Прикрепляю две реализации процесса. Файлы с подчёркиванием - im, без - re.

Простите, но в обоих представленных фрагментах SNR около 67дБ, SFDR чуть меньше 74дБ. Это действительно такие (явно не 100дб SNR) фрагменты?

Поделиться сообщением


Ссылка на сообщение
Поделиться на другие сайты

Вот, немного поправил программу, и подбавил аддитивного шума.

Теперь вторая синусоида в склейке разрешается при любых настройках, даже без коррекции её фазы.

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

% Формируем модель "реального" сигнала
signal = sin(2*pi*(0:4*len_buf)*freq/fsampl)+(1e-2)*randn(1, 4*len_buf+1)+...
     0.03*cos(2*pi*(0:4*len_buf)*(freq+0.5)/fsampl);

% Заполняем буферы (теперь всё точно)
buffer_1 = signal(1:len_buf);
sampl_gap = floor(1e6*gap/fsampl);
buffer_2 = signal(len_buf+1+sampl_gap:len_buf+sampl_gap+len_buf);

% Пытаемся склеить без учёта фазовых соотношений
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);

% Находим фазовые спектры кусков
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*len_buf*freq/fsampl, 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(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
title('Spectrum of short sequence')

figure (5)
semilogy (bins(1:1000)/2, abs(sp_sum_2(1:1000)))
ylim([0.1 1000]);
grid on
title('Spectrum of long sequence')
return;

post-4987-1217747456_thumb.jpg - модуль спектра последовательности одинарной длины,

 

post-4987-1217747466_thumb.jpg - модуль спектра склейки. B)

 

Продолжение следует. :)

Поделиться сообщением


Ссылка на сообщение
Поделиться на другие сайты

2 EKrishin

Скажите, в файлах, выложенных Вами, частота выходных отсчётов действительно 1 МГц? Или всё-таки менее?

Да, действительно, 1 Мгц. А что не так?

 

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

Да, именно это и интересно. Вообще хотелось получать выборку в, скажем, 65536 точек, строить большой спектр и смореть с разрешением порядка 15-20 Гц. Но аппаратные средства не позволяют. Поэтому ищем путь сделать нечто из различных кусков сигнала.

 

 

Простите, но в обоих представленных фрагментах SNR около 67дБ, SFDR чуть меньше 74дБ. Это действительно такие (явно не 100дб SNR) фрагменты?

Да, действительно. Указанные 100db - это, пожалуй, максимум. А данные две реализации - лишь пример.

Поделиться сообщением


Ссылка на сообщение
Поделиться на другие сайты

Фазировка блоков именно к этому и приводит - уровень паразитной подставки примерно -70дБ даже на модели.
На какой модели? Не могли бы обнародовать, а то непонятно, о чём это Вы?

 

Да, действительно, 1 Мгц. А что не так?
Нет, ничего...

 

...Вообще хотелось получать выборку в, скажем, 65536 точек, строить большой спектр и смореть с разрешением порядка 15-20 Гц. Но аппаратные средства не позволяют. Поэтому ищем путь сделать нечто из различных кусков сигнала.
Скажите, а почему предложение ув. alex_os не подходит (пост №71)?

Если не трудно, опишите систему сбора данных подробнее.

Поделиться сообщением


Ссылка на сообщение
Поделиться на другие сайты

Присоединяйтесь к обсуждению

Вы можете написать сейчас и зарегистрироваться позже. Если у вас есть аккаунт, авторизуйтесь, чтобы опубликовать от имени своего аккаунта.

Гость
К сожалению, ваш контент содержит запрещённые слова. Пожалуйста, отредактируйте контент, чтобы удалить выделенные ниже слова.
Ответить в этой теме...

×   Вставлено с форматированием.   Вставить как обычный текст

  Разрешено использовать не более 75 эмодзи.

×   Ваша ссылка была автоматически встроена.   Отображать как обычную ссылку

×   Ваш предыдущий контент был восстановлен.   Очистить редактор

×   Вы не можете вставлять изображения напрямую. Загружайте или вставляйте изображения по ссылке.

×
×
  • Создать...