ДПФ и дробные частоты.

[stalkernet]

Всем хорошего дня.

вобщем как говорится - инициатива всегда наказуема и хотелось как лучше - получилось как всегда...... посоветовал решение на свою голову.

Ситуация. Есть сигнал состоящий из каши частот от 1Гц до 50. длительность 1 сек.
Этот сигнал надо разложить по частотам от 1 до 50 с шагом 0,1 Гц. и получить реальную амплитуду.

поковырял инет. перепробывал различные предложения. ситуация не очень. кроме как 10 сек интервала измерения - достоверных данных не получил. все остальные предложения работают если сигнал с одной дробной частотой (типа 1, 1.5 и 2 Гц). Как только состав сигнала становится типа 1, 1.2, 1.5, 2 Гц все перестает правильно определять.

Теперь собственно вопрос( надеюсь грамотно сформулировал).

возьмем идеальную ситуацию.
Есть сигнал от 1 до 2Гц с шагом 0.1Гц. аплитуда и фаза одинаковы. длительность измерения 1 сек. частота оцифровки 10000 Гц. затухающих состовляющих нет.

требуется разложить на спект с шагом 0.1 Гц. и узнать реальную амплитуду каждой частоты. 1, 1.1, 1.2 и т.д. до 2 Гц.

Если есть идеи - буду рад. с обстенными пока TwoPi. тоесть 0.

[northener]

А при чём программирование?

[Pavia]

Короче. Сигналы бывают разных форм. Переводческие и непериодические.
Вообще любой непериодический сигнал мы можем разложить на переодический. В рамках допущей изложенных Фурье, о период гармонических сигналов должен быть одинаково кратным т.е T,2*T,3*T..n*T. Как струна зажатая с двух концов, которая принимает разные частоты.

Будем считать что в исходном сигнале периодические.
Если строить Фурье спектр. То он будет растекаться из-за того что во время записи мы оборвали его на скажем четверти периода. Когда как преобразование Фурье требует строгое совпадение периодов. Это уже не струна зажатая с двух сторон, а свободная.

Собственно с компенсировать это не проблема. Нам нужно всего лишь для каждой гармоники найти фазу. Идея такая же как и в алгоритме "микстура Гаусса"(Гауссова смесь распределений (GMM)). Перебираем все частоты, для частоты находим лучшую фазу. Путем расчёта корреляции. Где коэффициент корреляции выше тот и берем. После вычитаем найденную гармонику из исходного сигнала. Пока он не станет ниже некоторого порога.
И снова перебираем частоты и снова вычитаем пока неостанутся шумы.

А вообще 0,1 ГЦ замораживаться не стоит, а проще увеличить время измерения до 10 секунд. И сделать обычное БПФ.


Я фазу не перебирал. Но поиск частот делал вот таким кодом точность была хорошая до 0.001 Гц.

// Вычисления коэффициента Фурье для заданной частоты w
function FT(w:Real; N:NativeInt; a:PAReal):Complex;Overload;
var i:NativeInt;
N1:Real;
begin
Result.Re:=0;
Result.Im:=0;
N1:=1/N;
for i:=0 to N-1 do
 begin
 Result.Re:=Result.Re+A[i]*Cos(-2*Pi*i*w*N1);
 Result.Im:=Result.Im+A[i]*Sin(-2*Pi*i*w*N1);
 end;
Result.Re:=Result.Re*N1;
Result.Im:=Result.Im*N1;
end;


// Поиск частоты с максимальной амплитудой. В диапазоне f1 и f2
function DetectF(a:TArrayReal; dwSamplesPerSec:DWord; f1:Real=20; f2:Real=22500):Real;
var i:Integer;
 N:Integer;
 fd:Real;
 z:Complex;
 r,MaxR,f,MaxF:Real;
begin
n:=Length(a);
fd:=dwSamplesPerSec;
f:=Trunc(f1);
MaxR:=0;
repeat
z:=FT(f/(Fd/N),N,PAReal(@a[0]));
r:=Amp(z);
if r>MaxR then
  begin
  MaxF:=f;
  MaxR:=r;
  end;
f:=f+1;
until f>f2;

f:=MaxF-1;
MaxR:=0;
for i:=0 to 20 do
  begin
  z:=FT(f/(Fd/N),N,PAReal(@a[0]));
  r:=Amp(z);
  if r>MaxR then
    begin
    MaxF:=f;
    MaxR:=r;
    end;
  f:=f+0.1;
  end;

f:=MaxF-0.1;
MaxR:=0;
for i:=0 to 20 do
  begin
  z:=FT(f/(Fd/N),N,PAReal(@a[0]));
  r:=Amp(z);
  if r>MaxR then
    begin
    MaxF:=f;
    MaxR:=r;
    end;
  f:=f+0.01;
  end;

f:=MaxF-0.01;
MaxR:=0;
for i:=0 to 20 do
  begin
  z:=FT(f/(Fd/N),N,PAReal(@a[0]));
  r:=Amp(z);
  if r>MaxR then
    begin
    MaxF:=f;
    MaxR:=r;
    end;
  f:=f+0.001;
  end;

f:=MaxF-0.001;
MaxR:=0;
for i:=0 to 20 do
  begin
  z:=FT(f/(Fd/N),N,PAReal(@a[0]));
  r:=Amp(z);
  if r>MaxR then
    begin
    MaxF:=f;
    MaxR:=r;
    end;
  f:=f+0.0001;
  end;

result:=MaxF;
end;

// Поиск первой по частоте гармоники в сигнале. Выше некоторого уровня шума t. 
function DetectFirstF(a:TArrayReal; t:Real; dwSamplesPerSec:DWord; f1:Real=20; f2:Real=22500):Real;
var
 i,N,NewN:Integer;
 aa,b:TArrayReal;
 z:TArrayComplex;
 r:Real;
 fd:Real;
begin
fd:=dwSamplesPerSec;
N:=Length(a);
aa:=copy(a);
r:=Mean(aa);
sub(aa,r);        // Убираем постоянную составляющую.
CopyInRe(z,aa);

FFT(z,False);    // Получаем спектр частот

amp(b,z);
for i:=1 to Length(b)-2 do
 b[i]:=0.25*b[i-1]+0.5*b[i]+0.25*b[i+1]; // Убираем небольшие шумы, могут влиять на знак 1 производной.
for i:=6 to Length(b) div 2 do
 begin
 if (b[i]>t) and (b[i]-b[i-1]>0) and (b[i]-b[i+1]>0) then // проверяем порог и знак второй производной
   begin
   Result:=DetectF(aa,dwSamplesPerSec,i*(Fd/N)-5,i*(Fd/N)+5); // Уточняем частоту
   Exit;
   end;
 end;
end;

[stalkernet]

Pavia Огромное спасибо!!! всем бы так граммотно могли обьяснять.

Практически все понятно. думаю с адаптацией алгоритма вопросов не будет.
есть пару вопросов.

формат TArrayComplex - вещественное + мнимое или наобот (было и такое в алгоритмах)

CopyInRe(z,aa); - здесь копируем аа в вещественную часть?

на выходе процедуры FFT(z,False); получаем двух мерный массив (вещественное, мнимое), размерностью Fдискретизации /2 ? С дальнейшим приобразованием 2*Sqrt(Sqr(f[i].x) + Sqr(f[i].y))/NDiscretization

[Pavia]

Видно что Вы человек подготовленный. Поэтому с Вами проще.
Тут ка бы описано 2 разных алгоритма один текстом другой кодом. В коде там ищется первый локальный максимум, для реальных данных надо бы сгладить. Код писался 9 лет назад. Сейчас особо нет времени с этим разбираться возможно там частота дискретизации не учитывается и надо собственные частоты пересчитать на неё.

Complex=Packed record
Re,Im:Real;
end;
TArrayComplex=array of Complex;


CopyInRe(z,aa); - здесь копируем аа в вещественную часть массива, копируется один в другой.
FFT выдает как есть в начали положительные частоты потом отрицательные. С коэффициентом 1/N.
amp(b,z); уже делает Sqrt(Sqr(f[i].x) + Sqr(f[i].y))

Библиотечные функции можете можете тут взять
https://yadi.sk/d/y4Yzme0DPqKm4g

[Pavia]

...

[stalkernet]

Иван. в примере WaveAnaliz. нет файла dcvMap.pas DCU есть. но собрать под XE10. сам понимаешь
DCU только для определёной версии делфей.

если есть pas - скинь.

попытаюсь собрать через старые делфя типа V7 на виртуалке.

работаю на rad XE 10.4

во всяком случае не так крИтЫчно. ))) соберу на Д7. посмотрю результат.

Еще раз спасибо.