[Stanislav 0]

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

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

 

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

А вот и посмотрим, как она размажется. :)

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

 

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

А если хорошо подумать? ;)

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

 

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

Ну, это уже лишнее.

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

Это точно.

[alex_os 0]

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

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

[Stanislav 0]

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

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

 

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

 

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

 

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

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

2 EKrishin

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

Изменено 1 августа, 2008 пользователем Stanislav

[Stanislav 0]

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

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 раза! B)

Изменено 2 августа, 2008 пользователем Stanislav

[fontp 0]

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

 

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

 

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

Дайте две.

 

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

 

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

и подпись

 

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

 

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

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

[GetSmart 0]

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

и подпись

 

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

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

 

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

...

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

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

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

 

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

Изменено 2 августа, 2008 пользователем GetSmart

[Stanislav 0]

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

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

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

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

 

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

Дайте две.

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

 

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

 

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

и подпись

 

"барон Мюнхаузен" 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, наблюдайте внимательно:

 

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

 

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

 

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

 

 

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

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

Изменено 2 августа, 2008 пользователем Stanislav

[alex_os 0]

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

 

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

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

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

 

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

 

[Stanislav 0]

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

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

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

[alex_os 0]

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

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

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

[Stanislav 0]

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

Измерить!

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

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

[729 0]

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

 

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

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

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

[Stanislav 0]

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

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

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;

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

 

- модуль спектра склейки. B)

 

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

[EKirshin 0]

2 EKrishin

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

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

 

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

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

 

 

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

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

[Stanislav 0]

Фазировка блоков именно к этому и приводит - уровень паразитной подставки примерно -70дБ даже на модели.

На какой модели? Не могли бы обнародовать, а то непонятно, о чём это Вы?

 

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

Нет, ничего...

 

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

Скажите, а почему предложение ув. alex_os не подходит (пост №71)?

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