VPF::Цифровая фильтрация - Форум программистов

svip

Дата 27.3.2007, 16:49 (ссылка)

(нет голосов)

Загрузка ...

Шустрый

Профиль
Группа: Участник
Сообщений: 61
Регистрация: 27.3.2007

Репутация: нет
Всего: нет

Есть аналоговый сигнал частотой от 0 до 1 Гц. Оцифровали его через АЦП и нужно програмно выделить четыре частоты из него (0.12 0.8 0.6 0.4 Гц).
Помогите разобраться. (Курсовая горит)

svip

Дата 27.3.2007, 17:33 (ссылка)

(нет голосов)

Загрузка ...

Шустрый

Профиль
Группа: Участник
Сообщений: 61
Регистрация: 27.3.2007

Репутация: нет
Всего: нет

неужели никто не знает или нехочет помочь?

Alexeis

Дата 27.3.2007, 18:00 (ссылка)

(нет голосов)

Загрузка ...

Амеба

Профиль
Группа: Админ
Сообщений: 11743
Регистрация: 12.10.2005
Где: Зеленоград

Репутация: 109
Всего: 459

svip, легко прогони его через фурье (БПФ), в комплексном спектре оставь только нужную компоненту (слева и справа) и обратным БПФ получиться нужная составляющая.

Пример реализации фурье есть в DRKB.

Это сообщение отредактировал(а) Alexeis - 27.3.2007, 18:16

--------------------

Vit вечная память.

Обсуждение действий администрации форума производятся только в этом форуме

гениальность идеи состоит в том, что ее невозможно придумать

svip

Дата 28.3.2007, 20:25 (ссылка)

(нет голосов)

Загрузка ...

Шустрый

Профиль
Группа: Участник
Сообщений: 61
Регистрация: 27.3.2007

Репутация: нет
Всего: нет

а пример кода можно привести

Alexeis

Дата 28.3.2007, 22:19 (ссылка)

(нет голосов)

Загрузка ...

Амеба

Профиль
Группа: Админ
Сообщений: 11743
Регистрация: 12.10.2005
Где: Зеленоград

Репутация: 109
Всего: 459

Этот код не оптимальный так как использует вычисления с плавающей точкой и классической преобразование фурье для 4096 точек, частота дискретизации принята 400 Гц.

Внимание, расчет достаточно длительный и на процессоре Athlon 3800+ занял несколько секунд.

Для столь низких частот желательно уменьшить частоту дискретизации например 40Гц, и соответственно уменьшить число точек преобразования например до 512, тогда скорость многократно возрастет. Длинные буфера нужно разбивать на блоки по длине равной длине буфера фурье и преобразовывать поблочно.

Вот такой получился простенький пример

Код


const
  Diskr  = 400; //выборок в сек
  fr1 = 0.12;
  fr2 = 0.6;
  fr3 = 0.8;
  fr4 = 0.4;

procedure ClassicDirect(var aSignal, aSpR, aSpI: array of Double; N: LongInt);
var
  lSrch : LongInt;
  lGarm : LongInt;
  dSumR : Double;
  dSumI : Double;
begin
  for lGarm := 0 to N div 2 - 1
  do
     begin
       dSumR := 0;
       dSumI := 0;

       for lSrch := 0 to N - 1
       do
          begin
            dSumR := dSumR + aSignal[lSrch] * Cos(lGarm * lSrch / N * 2 * PI);
            dSumI := dSumI + aSignal[lSrch] * Sin(lGarm * lSrch / N * 2 * PI);
          end;

       aSpR[lGarm] := dSumR;
       aSpI[lGarm] := dSumI;
     end;
end;

procedure ClassicInverce(var aSpR, aSpI, aSignal: array of Double; N: LongInt);
var
  lSrch : LongInt;
  lGarm : LongInt;
  dSum  : Double;

begin
  for lSrch := 0 to N - 1
  do
     begin
       dSum := 0;
       for lGarm := 0 to N div 2 - 1
       do
         dSum := dSum + aSpR[lGarm] * Cos(lSrch * lGarm * 2 * Pi / N)
                      + aSpI[lGarm] * Sin(lSrch * lGarm * 2 * Pi / N);
       aSignal[lSrch] := dSum * 2 / N;
     end;
end;


function Signal(n : integer) : double;
begin
  Result := 20 * Sin (fr1 * n / Diskr * 2 * Pi) +
            17 * Sin (fr2 * n / Diskr * 2 * Pi) +
            18 * Sin (fr3 * n / Diskr * 2 * Pi) +
            14 * Sin (fr4 * n / Diskr * 2 * Pi);
end;

{$O-}

procedure TForm1.Button1Click(Sender: TObject);
var
  Source : array[0..4095] of Double;
  Cx, Cy : array[0..4095] of Double;
  Fr0_4  : array[0..4095] of Double;
   Gr   : array[0..4095] of Double;
  i, n : integer;

begin
  Image1.Canvas.MoveTo(i, Round(100 - Source[0]));
  for I := 0 to 4095
  do
    Begin
      Source[i] := Signal(i);
      if (i mod 5) = 0
      then
        Image1.Canvas.LineTo(Round(i / 5), Round(100 - Source[i]));
    End;

  ClassicDirect(Source, Cx, Cy, 4096);
  n := round(fr1 / Diskr * 4096);
 {    // спектр
  Image1.Canvas.MoveTo(0, Round(300 - Sqrt(Sqr(Cx[0]) + Sqr(Cy[0])) / 200));

  for I := 0 to 40
  do
    Begin
      Gr[i] := Sqrt(Sqr(Cx[i]) + Sqr(Cy[i])) / 200;
      Image2.Canvas.LineTo(i * 10, Round(300 - Gr[i]));
    End;
    }

  for I := 0 to 2047
  do
    if abs(i - n) > 1
    then
      Begin
        Cx[i] := 0;
        Cy[i] := 0;
      End;
  ClassicInverce(Cx, Cy, Fr0_4, 4096);

  Image2.Canvas.MoveTo(0, Round(100 - Fr0_4[0]));
  for I := 1 to 4095
  do
    if (i mod 5) = 0
    then
      Image2.Canvas.LineTo(Round(i / 5), Round(100 - Fr0_4[i]));

  i := i;
end;

--------------------

Alexeis

Дата 29.3.2007, 09:44 (ссылка)

(нет голосов)

Загрузка ...

Амеба

Профиль
Группа: Админ
Сообщений: 11743
Регистрация: 12.10.2005
Где: Зеленоград

Репутация: 109
Всего: 459

Вот наконец нашел реализацию алгоритма БПФ для целочисленных операций (к сожалению нет сведений об авторе кода)

Алгоритм быстрого преобразование Фурье на 1024 точек (перевод в частотную область)
//вспомогательные массивы для ускорения обработки

Код


var int1024sin2pidiv1024, int1024cos2pidiv1024: array[0..1023] of integer; 
for i:=0 to 1023 do begin  //массивы заполняются один раз при инициализации.
int1024cos2pidiv1024[i]:=round(1024*cos(2*pi*i/1024)); //один период косинуса:1024 отсчетов        int1024sin2pidiv1024[i]:=round(1024*sin(2*pi*i/1024)); //один период синуса:1024 отсчетов        
end;
var re,im : array[0..1023] of integer; //re[], im[] – действительная и мнимая части сигнала
// re[] заполняется отсчетами сигнала (в обычном порядке) до вызова процедуры
// im[] заполняется нулями до вызова процедуры

procedure FFT1024();
var i,j,k,L,N, ind, md2,mdd,mdd1, i1,i2:integer;
    a,b,c,d,batr,bati:integer;
    nm:array[0..1023] of integer;
begin
  N:=1024;
  L:=10;         //L:=round(ln(N)/ln(2));
  ind:=0;        //перестановка отсчетов в битреверсном порядке как для re[] так и для im[]
  nm[0]:=0;
  for i:=0 to L-1 do begin
    for j:=0 to (1 shl i)-1 do begin
      if ind>nm[ind] then begin
        a:=re[ind];                                \\
        re[ind]:=re[nm[ind]];               \\ перестановка для re[]
        re[nm[ind]]:=a;                        \\
        a:=im[ind];                               \\ 
        im[ind]:=im[nm[ind]];             \\ перестановка для im[]
        im[nm[ind]]:=a;                       \\
      end;
      ind:=ind+1;
      nm[ind]:=nm[j]+1 shl (L-1-i);
    end;
  end;     //конец перестановки отсчетов в битреверсном порядке
  md2:=N shr 1;
  for k:=0 to L-1 do begin     //цикл по стадиям  (10 стадий)
    mdd:=md2 shr k;
    mdd1:=md2 shr(L-k-1);
    for i:=0 to mdd-1 do begin  //цикл по группам в каждой стадии  (512x1, 256x2, 128x4, 64x8, 16x32,…, 1x512)
      i1:=2*mdd1*i;
      i2:=2*mdd1*i+mdd1;
      for j:=0 to mdd1-1 do begin   //цикл по парам внутри группы (выполнение бабочки для каждой пары в группе)
        a:=re[j+i1];
        b:=im[j+i1];
        c:=re[j+i2];
        d:=im[j+i2];
        ind:=j*mdd;
        batr:=(c*int1024cos2pidiv1024[ind] + d*int1024sin2pidiv1024[ind])div 1024; 
        bati:=(d*int1024cos2pidiv1024[ind] - c*int1024sin2pidiv1024[ind])div 1024;  
        re[j+i1]:=a+batr;     \\бабочка
        im[j+i1]:=b+bati;    \\бабочка
        re[j+i2]:=a-batr;      \\бабочка
        im[j+i2]:=b-bati;     \\бабочка
      end;
    end;
  end;
  for k:=0 to N-1 do begin  //деление на 1024 (арифметический сдвиг вправо на 10 дв. разрядов)
    re[k]:=re[k] div N;
    im[k]:=im[k] div N;
  end;
end;

Код


//Алгоритм быстрого обратного преобразования Фурье на 1024 точек (перевод спектра в временнýю область)

procedure IFFT1024();
var i,j,k,L,N, ind, md2,mdd,mdd1,
    i1,i2:integer;
    a,b,c,d,batr,bati:integer;
    nm:array[0..1023] of integer;

begin
  N:=1024;
  L:=10; //L:=round(ln(N)/ln(2));

  ind:=0;            //перестановка отсчетов в битреверсном порядке как для re[] так и для im[]
  nm[0]:=0;
  for i:=0 to L-1 do begin                       
    for j:=0 to (1 shl i)-1 do begin
      if ind>nm[ind] then begin
        a:=re[ind];
        re[ind]:=re[nm[ind]];
        re[nm[ind]]:=a;
        a:=im[ind];
        im[ind]:=im[nm[ind]];
        im[nm[ind]]:=a;
      end;
      ind:=ind+1;
      nm[ind]:=nm[j]+1 shl (L-1-i);
    end;
  end;                 //конец перестановки отсчетов в битреверсном порядке

  md2:=N shr 1;
  for k:=0 to L-1 do begin         //цикл по стадиям  (10 стадий)
    mdd:=md2 shr k;
    mdd1:=md2 shr(L-k-1);
    for i:=0 to mdd-1 do begin   //цикл по группам в каждой стадии  (512x1, 256x2, 128x4, 64x8, 16x32,…, 1x512)    
      i1:=2*mdd1*i;
      i2:=2*mdd1*i+mdd1;
      for j:=0 to mdd1-1 do begin  //цикл по парам внутри группы (выполнение бабочки для каждой пары в группе)
        a:=re[j+i1];
        b:=im[j+i1];
        c:=re[j+i2];
        d:=im[j+i2];

  ind:=j*mdd;
  batr:=(c*int1024cos2pidiv1024[ind] - d*int1024sin2pidiv1024[ind])div 1024;
  bati:=(d*int1024cos2pidiv1024[ind] + c*int1024sin2pidiv1024[ind])div 1024;

        re[j+i1]:=a+batr;
        im[j+i1]:=b+bati;
        re[j+i2]:=a-batr;
        im[j+i2]:=b-bati;
      end;
    end;
  end;
end;

--------------------

svip

Дата 29.3.2007, 20:43 (ссылка)

(нет голосов)

Загрузка ...

Шустрый

Профиль
Группа: Участник
Сообщений: 61
Регистрация: 27.3.2007

Репутация: нет
Всего: нет

ограмное спасибо за примеры.
На сколько я понял
  fr1 = 0.12;
  fr2 = 0.6;
  fr3 = 0.8;
  fr4 = 0.4;

мы складываем и получаем смешанную частоту, а что получаем путем вычисления?

со вторым примеров вообще не разобрался

если можете объясните поподробнее на пальцах непонятливому студенту

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

Alexeis

Дата 29.3.2007, 23:24 (ссылка)

(нет голосов)

Загрузка ...

Амеба

Профиль
Группа: Админ
Сообщений: 11743
Регистрация: 12.10.2005
Где: Зеленоград

Репутация: 109
Всего: 459

Цитата(svip @ 29.3.2007, 20:43

)

мы складываем и получаем смешанную частоту,

У меня же нет оцифрованного сигнала. Вот я его заменил искуственным. Путем вычисления сначала получаем комплексный спектр сигнала, затем оставляем все компоненты кроме нужной нам, а остальные зануляем. После чего обратным преобразованием фурье получем отфльтрованную компоненту (то что требуется по условию). Но приведенный алгоритм классического фурье медленный, потому второй пример есть реализация целочисленного БПФ, который многократно быстрее, но с ним и больше возни.

Код


 Image1.Canvas.MoveTo(i, Round(100 - Source[0]));
  for I := 0 to 4095
  do
    Begin
      Source[i] := Signal(i);
      if (i mod 5) = 0
      then
        Image1.Canvas.LineTo(Round(i / 5), Round(100 - Source[i]));
    End;

Рисуем исходный сигнал.

ClassicDirect(Source, Cx, Cy, 4096); - прямое фурье
n := round(fr1 / Diskr * 4096); - номер нужной компоненты в спектре

Код


for I := 0 to 2047
  do
    if abs(i - n) > 1
    then
      Begin
        Cx[i] := 0;
        Cy[i] := 0;
      End;

Удаление лишних компонент

ClassicInverce(Cx, Cy, Fr0_4, 4096); восстановление сигнала (обратное фурье)

Код


 for I := 1 to 4095
  do
    if (i mod 5) = 0
    then
      Image2.Canvas.LineTo(Round(i / 5), Round(100 - Fr0_4[i]));

Рисуем результат

--------------------

svip

Дата 30.3.2007, 15:42 (ссылка)

(нет голосов)

Загрузка ...

Шустрый

Профиль
Группа: Участник
Сообщений: 61
Регистрация: 27.3.2007

Репутация: нет
Всего: нет

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

Alexeis

Дата 30.3.2007, 16:19 (ссылка)

(нет голосов)

Загрузка ...

Амеба

Профиль
Группа: Админ
Сообщений: 11743
Регистрация: 12.10.2005
Где: Зеленоград

Репутация: 109
Всего: 459

Точно также как и обычное фурье, только для его работы сначала вычисляется таблица синусов и косинусов и данные должны быть целочисленными.
re,im - входные выходные массивы данных как для прямого так и для обратного фурье.
Для прямого исходные данные помещаются в re, im должен быть заполнен нулями.
Для обратного входные это re, im, результат возвращается в re. Если что его можно "подрихтовать" для нужного количества точек.

--------------------

0 Пользователей читают эту тему (0 Гостей и 0 Скрытых Пользователей)
0 Пользователей:
« Предыдущая тема \| Delphi: Общие вопросы \| Следующая тема »