Проблемы с реализацией алгоритма Герцеля

IONEX
Сообщения: 5
Зарегистрирован: 02 ноя 2011, 00:42

Проблемы с реализацией алгоритма Герцеля

Сообщение IONEX »

Доброго всем времени суток!

Задача: имею вектор данных (реальная измерительная информация) допустим 3000 отсчетов. Темп данных 1 Гц. Из всей смеси шумов и сигнала необходимо выделить только сигнал с периодом 400 сек. Как подсказали на форуме для этого необходимо использование алгоритма Герцеля. Программу реализовал на MATLAB по примеру кода размещенного на сайте написанного в коде С (тестовые примеры сошлись, поэтому здесь ошибки нет).

Проблема: Хочу выполнить проверку. Для этого создал синусоиду длительностью 3000 сек (темп данных 1сек или 1Гц) с периодом 400 секунд.

[img]
http://files.mail.ru/95F4IY?t=1
[/img]

Ожидаю на выходе значение амплитуды 1 для k=400, N=3000. Для этого в конце процедуры поставил расчет модуля

Код: Выделить всё

%НА MATLABе
level=abs(SR-1i*SI);
Ответ выдает level=2.46, для других периодов k=100, 200 и т.д. ожидаю получение level=0. (машинный ноль), однако при этом значения level далеки от 0 и могут быть равными 26.

Подскажите в чем может быть ошибка?

ivan219
Сообщения: 61
Зарегистрирован: 09 май 2011, 16:39

Re: Проблемы с реализацией алгоритма Герцеля

Сообщение ivan219 »

K - это частота искомого сигнала
N- длинна иследуемого сигнала.
Если вы хотите искать значение только одно частоты во всей выборке то k не должно меняться.

IONEX
Сообщения: 5
Зарегистрирован: 02 ноя 2011, 00:42

Re: Проблемы с реализацией алгоритма Герцеля

Сообщение IONEX »

Ну если мне нужно найти гормоники на двух частототах, почему я не могу на вход процедуры подавать необходимое мне к (например найти гармонику с частотой 1/400 и 1/200 из данных исходного сигнала? При этом каждый раз буду все расчитывать заново с нужным мне к)?

Вопрос состоял в том: "Почему не получаються ожидыемые мной значения для тестового гармонического сигнала с периодом 400 сек"

1) k=400, то при расчете модуля полученных комплексных коэффициентов амплитуда не равна 1.
2) k=200, казалось бы, что поскольку одна полоска в спектре на 1/400 Гц, то для k=200 модуль будет равен 0, а он может достигать амплитуды 26 (для других к, отличных от 400 значения могут быть выше).

Привожу интерпритированный мной код с С на MATLAB. В коде представленом ниже добавлены только строки которые отслеживаю пропуски в данных (пропуск соотвествует данных == 0)

Код: Выделить всё

function [level]=hertzel_filter(data,k)
%Процедура выделения гармоники k из массива данных data с использованием 
%алгоритма Герцеля.
%Входные данные:
%    data - массив данных для которых необходимо выделить гармонику
%    к - частота гармоники из учета что data это герцовки
%Выходные данные:
%    level - амплитуда искомой гармоники

[N,y]=size(data);

alpha=2*cos(2*pi*k/N);

%Расчет поворотных коэффициернтов (мнимая и реальные части)
wr=cos(2*pi*k/N); %реальная часть
wi=sin(2*pi*k/N); %мнимая часть

iteration=zeros(N,y);%Матрица промежуточных результатов
for sv=1:y
    timeOBS=find(data(:,sv));%Поиск данных
    if ~isempty(timeOBS)%Если данные есть, то:
        indexD=diff(timeOBS);%поиск пропусков в данных
        temp=find(indexD>300);%Если разрыв между данными более 300 сек, это новый интервал (определение пропусков выше 300 сек.)
        %--------------Определение начала и конца измерений----------------
        pf=zeros(1,1+length(temp));
        zf=pf;
        pf(1)=timeOBS(1);
        pf(2:end)=timeOBS(temp+1);
        zf(1:end-1)=timeOBS(temp);
        zf(end)=timeOBS(end);
        for m=1:1:1+length(temp)
            %---------Учет первых элементов промежуточной матрицы----------
            iteration(pf(m),sv)=data(pf(m),sv);
            iteration(pf(m)+1,sv)=data(pf(m)+1,sv)+alpha*data(pf(m),sv);
            %------Итерационный процесс расчета промежуточной матрицы------
            for i=pf(m)+2:zf(m)
                iteration(i,sv)=data(i,sv)+alpha*iteration(i-1,sv)-iteration(i-2,sv);
            end
            SR=iteration(zf(m),sv)*wr-iteration(zf(m)-1,sv);
            SI=iteration(zf(m),sv)*wi;
            level=abs(SR-1i*SI);
        end
    end
end

ivan219
Сообщения: 61
Зарегистрирован: 09 май 2011, 16:39

Re: Проблемы с реализацией алгоритма Герцеля

Сообщение ivan219 »

Код: Выделить всё

level=abs(SR-1i*SI);
Возможно ошибаюсь не знаю синтаксис Матлаба но модуль рассчитывается так

Код: Выделить всё

level=Sqrt(Sqr(SR) + Sqr(SI)) / N;

Аватара пользователя
Бахурин Сергей
Администратор
Сообщения: 1114
Зарегистрирован: 05 окт 2010, 19:55
Контактная информация:

Re: Проблемы с реализацией алгоритма Герцеля

Сообщение Бахурин Сергей »

IONEX писал(а):Подскажите в чем может быть ошибка?

Вот код для матлаба для алгоритма Герцеля думаю вам поможет

Код: Выделить всё


Fs = 1;
f0 = 1/400;
N = 3000;
t = (0:N-1)/Fs;

k = fix(N*f0/Fs);


s = cos(2*pi*f0*t);
alpha = 2*cos(2*pi*k/N);

W = exp(i*2*pi*k/N);

a = [1 -alpha 1];

S = filter(1,a,s);

S0 = fft(s);
SQ = abs(S0(k+1))
SR = abs(S(N)*W-S(N-1))

результаты SQ и SR полностью совпадают

IONEX
Сообщения: 5
Зарегистрирован: 02 ноя 2011, 00:42

Re: Проблемы с реализацией алгоритма Герцеля

Сообщение IONEX »

Хорошо БПФ и Герцель сошлись, это доказывает правильность работы. Непойму заначение результата который они вывели (около 1000), чему соответсвует полученное знание (амплитуда гармоники на частоте f0 ведь равна 1)?

Аватара пользователя
Бахурин Сергей
Администратор
Сообщения: 1114
Зарегистрирован: 05 окт 2010, 19:55
Контактная информация:

Re: Проблемы с реализацией алгоритма Герцеля

Сообщение Бахурин Сергей »

IONEX писал(а):Хорошо БПФ и Герцель сошлись, это доказывает правильность работы. Непойму заначение результата который они вывели (около 1000), чему соответсвует полученное знание (амплитуда гармоники на частоте f0 ведь равна 1)?
Если говорить строго то амплитуда спектра синусоиды равна 0.5 на частоте синусоиды если амплитуда синусоиды 1. ОДнако при дискретном преобразовании Фурье вы считаете амплитуду не в вольтах а в абстрактных числах. Можете сделать нормировку, т.е. поделить на N и тогда получите примерно 0.333. Почему не 0.5? Потому что вы произвели расчет на частоте не попадающей на бин спектра и ваш спектр растекся. Про этот эффект вы можете почитать здесь

IONEX
Сообщения: 5
Зарегистрирован: 02 ноя 2011, 00:42

Re: Проблемы с реализацией алгоритма Герцеля

Сообщение IONEX »

Да, согласен, 0.5. Написал, только потом заметил. Спасибо за помощь, кажется, это то что мне было необходимо.


По мере написания кода возник еще один вопрос: как предусмотреть пропуск данных (в векторе данных содержится 0), весовых матриц нет, как например это есть при МНК оценке?

Аватара пользователя
Бахурин Сергей
Администратор
Сообщения: 1114
Зарегистрирован: 05 окт 2010, 19:55
Контактная информация:

Re: Проблемы с реализацией алгоритма Герцеля

Сообщение Бахурин Сергей »

IONEX писал(а): По мере написания кода возник еще один вопрос: как предусмотреть пропуск данных (в векторе данных содержится 0), весовых матриц нет, как например это есть при МНК оценке?
что вы имеете ввиду под пропуском данных?

IONEX
Сообщения: 5
Зарегистрирован: 02 ноя 2011, 00:42

Re: Проблемы с реализацией алгоритма Герцеля

Сообщение IONEX »

Это реальная измерительная информация (навигационные сообщения) получемая с бортов спутников системы GPS. Есть такие моменты времени, когда слежение за спутником пропадает (затенение спутника зданиями и др. объектами), при этом таким данным присваивается значение 0, а затем снова возобновляется и приемник продолжает принимать сигналы дальше.
Таким образом, я имею данные в которых могут содержаться пропуски наблюдений (ноль в векторе данных). Теперь мне из всего вектора данных (информация + шум) необходимо выделить амплитуду гармоники с периодом 400 сек (что описано выше). Вопрос заключается в том, что поскольку фильтр использует обратную связь, а в векторе данных могут существовать пропуски, это наверника приведет к неправильному результату, как можно исключить влияние пропусков на результат?

Пока писал ответ возникло пару идей как это преодолеть:
1 Это использование фильтра без ОС (реализацию таких фильтров пока не знаю)
2 Поскольку пропуски данных сравнительно небольшие (десятки..сотни секунд) по отношению к интервалу видимости спутника (десятки минут.. шесть часов) возможна простая "склейка" (удаление нулей из их позиций и смещение всех данных влево), а затем производить фильтрацию с использованием алгоритма Герцеля.

Ответить