Ну если мне нужно найти гормоники на двух частототах, почему я не могу на вход процедуры подавать необходимое мне к (например найти гармонику с частотой 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