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

svip

Дата 13.10.2008, 20:04 (ссылка)

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

Загрузка ...

Шустрый

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

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

Есть сигнал 20Гц - 50 КГц

из него нужно отфильтровать полосу шириной 5 Гц (то есть 25-30Гц или 101-106Гц и тд.)

Значения амплитуды сигнала лежат в массиве.
Читал что нужно выполнить прямое преобразование Фурье, потом что-то как-то сделать, потом выполнить обратное преобразование фурье.
Но вот что-то ничего толком непонятно. Может у кого есть пример фильтра полосового??? или кто-то может понятно "на пальцах" с примерами объяснить. Буду очень благодарен, потому что уже пропло чувство, что наконец то смогу сделать этот фильтр.

svip

Дата 16.10.2008, 18:35 (ссылка)

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

Загрузка ...

Шустрый

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

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

народ помогите.

Alexeis

Дата 16.10.2008, 19:27 (ссылка)

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

Загрузка ...

Амеба

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

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

svip, прямой фурье дает комплексный спектр. т.е. комплексные амплитуды частот. Из спектра удаляются все частоты, т.е. зануляются соответствующие элементы массивов Re и Im, потом над Re и Im массивами производиться обратное преобразование фурье, вот и вся задача. У меня есть готовое решение, но каждый раз забываю. Стукни днем в асю, напомни smile

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

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

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

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

Alexeis

Дата 17.10.2008, 10:15 (ссылка)

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

Загрузка ...

Амеба

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

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

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

Также выкладываю целочисленный БПФ на 1024 точки "Бабочка".

Код


unit Unit1;

interface

uses
  Windows, Messages, SysUtils, Variants, Classes, Graphics, Controls, Forms,
  Dialogs, ExtCtrls, StdCtrls;

type
  TForm1 = class(TForm)
    btn1: TButton;
    Memo1: TMemo;
    imgpb1: TImage;
    procedure FormCreate(Sender: TObject);
    procedure btn1Click(Sender: TObject);
  private
    { Private declarations }
  public
    { Public declarations }
  end;

var
  Form1: TForm1;

implementation

{$R *.dfm}

uses math;

const
   MaxSize = 512;
   Ampl = 511;

var
 int1024sin2pidiv1024,
 int1024cos2pidiv1024: array[0..MaxSize - 1] of integer;
 re,
 im                  : array[0..MaxSize - 1] of integer;

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..MaxSize - 1] of integer;
begin
  N:=MaxSize;
  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 Ampl;
        bati:=(d*int1024cos2pidiv1024[ind] -
               c*int1024sin2pidiv1024[ind])div Ampl;
        re[j+i1]:=a+batr;     //бабочка
        im[j+i1]:=b+bati;    //бабочка
        re[j+i2]:=a-batr;      //бабочка
        im[j+i2]:=b-bati;     //бабочка
      end;
    end;
  end;

end;


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..MaxSize - 1] of integer;

begin
  N:= MaxSize;
  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 Ampl;
        bati:=(d*int1024cos2pidiv1024[ind] +
               c*int1024sin2pidiv1024[ind]) div Ampl;

        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;

   for i := 0 to MaxSize - 1
    do
      form1.memo1.Lines.Add(IntTostr(re[i]));
    form1.memo1.Lines.SaveToFile('c:\delpi.txt');

end;


procedure TForm1.FormCreate(Sender: TObject);
var
  i : integer;
begin
  for i:=0 to MaxSize - 1
  do
    begin  //массивы заполняются один раз при инициализации.
      int1024cos2pidiv1024[i]:=round(Ampl*cos(2*pi*i/MaxSize)); //один период косинуса:1024 отсчетов
      int1024sin2pidiv1024[i]:=round(Ampl*sin(2*pi*i/MaxSize)); //один период синуса:1024 отсчетов
      re[i] := round(100 * exp(-power(i - 200, 2) / 80) *
                           sin(40 * 2 * pi * i / MaxSize) -
                      80 * exp(-power(i - 210, 2) / 20) *
                           sin(45 * 2 * pi * i / MaxSize) -
                      30 * exp(-power(i - 205, 2) / 10) *
                           sin(41.5 * 2 * pi * i / MaxSize));
      im[i] := 0;
    end;
end;



procedure TForm1.btn1Click(Sender: TObject);
var
  i : integer;
begin
  imgpb1.Canvas.Pen.Color := $0000FF;
  imgpb1.Canvas.MoveTo(0, 120);
  for i := 0 to MaxSize - 1
  do
    imgpb1.Canvas.LineTo(i, 120 - re[i]);

  FFT1024();

{  for i := 0 to MaxSize - 1
  do
    begin
      re[i] := 0;
      im[i] := 0;
    end;
                  }
  //re[4] := 10000;

  IFFT1024();

  imgpb1.Canvas.Pen.Color := $009900;
  imgpb1.Canvas.MoveTo(0, 350);
  for i := 0 to MaxSize - 1
  do
    imgpb1.Canvas.LineTo(i, 400 {350} - re[i]);
                  // round(sqrt(re[i] * re[i] + im[i]*im[i])*2) );
  {
    for i := 0 to MaxSize - 1
    do
      memo1.Lines.Add(IntTostr(re[i]) + '  ' +
                      intToStr(im[i]));
    memo1.Lines.SaveToFile('c:\delpi.txt');
    }
end;

end.

Если что есть и рекурсивный фильтр, он быстрее работает.

Присоединённый файл ( Кол-во скачиваний: 59 )

____________.rar 184,78 Kb

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

svip

Дата 21.10.2008, 12:59 (ссылка)

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

Загрузка ...

Шустрый

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

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

Отлично то что нужно. буду разбираться, но сразу два вопроса, что за время меряется при клике на кнопку test, и в каких частотных пределах все это работает.

Alexeis

Дата 21.10.2008, 14:07 (ссылка)

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

Загрузка ...

Амеба

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

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

Цитата(svip @ 21.10.2008, 11:59

)

Ну там я сравнивал целочисленное фурье и с плавающей точнкой, за много итераций. Чем меньше время тем быстрее выполняется.
Частота полос среза зависит от частоты выборок. Тут просто массив, если его воспроизводить с разной скоростью, то частота будет меняться. Самая крайняя комплексная амплитуда соответствует половине частоты дискретизации, крайняя слева частоте волны, которая укладывается целиком в интервал над которым производиться фурье, т.е. будет зависеть от количества точек фурье.

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

svip

Дата 23.10.2008, 21:16 (ссылка)

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

Загрузка ...

Шустрый

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

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

Тут просто массив, если его воспроизводить с разной скоростью, то частота будет меняться.
А как воспроизводить с разной скоростью??

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

А тут вообще не понял.

svip

Дата 23.10.2008, 22:07 (ссылка)

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

Загрузка ...

Шустрый

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

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

И вот тоже не могу разобраться
40 * Sin(2 * pi / 100 * x)

тут 40- амплитуда,
x- время (допустим в секундах)
1/100 - частота, только вот в чем???

и как определить частоту фильтрации, получается есть использовать формулу 40 * Sin(2 * pi / 10 * x)
то пик амплитуды АЧХ лежит между SpeedButton1.left=2 и SpeedButton2.left=3

каксвязать с частотой??? чтобы у АЧХ проградуировать шкалу частот.

Alexeis

Дата 23.10.2008, 23:08 (ссылка)

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

Загрузка ...

Амеба

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

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

Цитата(svip @ 23.10.2008, 20:16

)

Тут просто массив, если его воспроизводить с разной скоростью, то частота будет меняться.
А как воспроизводить с разной скоростью??

Ну как, задать частоту выборок звуковой карте, она соответственно будет ее выдерживать. например 8000 выборок в секунду или 44100, тон будет разным.

Цитата(svip @ 23.10.2008, 20:16

)

А тут вообще не понял.

Спектр видел? Это ж массив значений. Рабочим является первая половина массива. Нулевой элемент амплитуда нулевой частоты, т.е. постоянная составляющая. Частота первого элемента спектра это величина обратная длительности всего участка (например 1024 точек) т.е. 1/T. Второго элемента 2/T, третьего 3/T т.д. В середине массива длительность 2х выборок, т.е. половина частоты дискретизации. Для 44100 это будет 22050Гц.
Советую немного почитать по дискретному фурье, мож яснее станет.

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

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