unit FIRLPF;

interface

uses
  Windows, Sysutils, Classes, Math;

type
  TDoubleArray = Array of Double;
  PDoubleArray = ^TDoubleArray;
  TWindowFunc = procedure (Count:Cardinal; Buffer:PDoubleArray);
  TProgressFunc = procedure (Progress, Max : Integer); register;

const
  Win_Kaiser = 0;
  Win_Blackman = 1;
  Win_Hamming = 2;
  Win_Hanning = 3;
  Win_Sin = 4;
  Win_Rectangler = 5;

procedure GetFilterCoef(Count, SamplingFreq, CutOffFreq:Cardinal; Buffer:PDoubleArray; WindowType:Word);
function FilteringProcess(Filter:PDoubleArray; FilterSize:Cardinal; InData:TList):Integer;
{procedure LowPassFilterForWaveFile(Src:String; Dest:String; CutOffFreq:Cardinal;
                                          WindowType:Word; ProgressFunc:TProgressFunc = nil);
}
procedure Window_Kaiser(Count:Cardinal; Buffer:PDoubleArray);
procedure Window_Hanning(Count:Cardinal; Buffer:PDoubleArray);
procedure Window_Hamming(Count:Cardinal; Buffer:PDoubleArray);
procedure Window_Blackman(Count:Cardinal; Buffer:PDoubleArray);
procedure Window_Sin(Count:Cardinal; Buffer:PDoubleArray);
procedure Window_Rectangler(Count:Cardinal; Buffer:PDoubleArray);

implementation

procedure Window_Kaiser(Count:Cardinal; Buffer:PDoubleArray);

  // 減衰率算出
  function Alpha(A:Byte):Extended;
  begin
      Result := 0;

      case A of
       0..21   : Result := 0;
       22..49  : Result := 0.5842*Power((A-21), 0.4) + 0.07886 * (A-21);
       50..255 : Result := 0.1102*(A-8.7);
      end;
  end;

  // 階乗
  function Factorial(n:Word):Cardinal;
  var
      i : Integer;
  begin
      Result := 1;

      for i := n downto 2 do
      begin
            Result := Result * i;
      end;
  end;

  // 0次第1種変形ベッセル関数
  function Bessel(u:Extended):Extended;
  const
      LOOP_LIMIT = 17;
  var
      i : Integer;
  begin
      Result := 1;

      for i := 1 to LOOP_LIMIT do
      begin
            Result := Result + Power( Power((u/2),i) / Factorial(i) , 2);
      end;
  end;

const
      Attenuation = 90;
var
      i : Integer;
      a : Double;
      w : Extended;
begin
      a := Alpha(Attenuation);

      for i := 0 to Count-1 do
      begin
            w := Bessel(A*sqrt(1-Power(i/(Count-1),2))) / Bessel(A);
            PDouble(Buffer)^ := PDouble(Buffer)^ * w;
            Inc(PDouble(Buffer));
      end;
end;

procedure Window_Hanning(Count:Cardinal; Buffer:PDoubleArray);
var
      i : Integer;
      w : Double;
begin
      for i := 0 to Count-1 do
      begin
            w := (0.5-0.5*cos((2*pi*i)/(Count-1)));
            PDouble(Buffer)^ := PDouble(Buffer)^ * w;
            Inc(PDouble(Buffer));
      end;
end;

procedure Window_Hamming(Count:Cardinal; Buffer:PDoubleArray);
var
      i : Integer;
      w : Double;
begin
      for i := 0 to Count-1 do
      begin
            w := (25/46-(1-25/46)*cos((2*pi*i)/(Count-1)));
            PDouble(Buffer)^ := PDouble(Buffer)^ * w;
            Inc(PDouble(Buffer));
      end;
end;

procedure Window_Blackman(Count:Cardinal; Buffer:PDoubleArray);
var
      i : Integer;
      w : Double;
begin
      for i := 0 to Count-1 do
      begin
            w := 0.423-0.498*cos((2*pi*i)/(Count-1))+0.0792*cos((4*pi*i)/(Count-1));
            PDouble(Buffer)^ := PDouble(Buffer)^ * w;
            Inc(PDouble(Buffer));
      end;
end;

procedure Window_Sin(Count:Cardinal; Buffer:PDoubleArray);
var
      i : Integer;
      w : Double;
begin
      for i := 0 to Count-1 do
      begin
            w := sin(pi*i/(Count-1));
            PDouble(Buffer)^ := PDouble(Buffer)^ * w;
            Inc(PDouble(Buffer));
      end;
end;

procedure Window_Rectangler(Count:Cardinal; Buffer:PDoubleArray);
begin
      // レクタンギュラは何もなっしんぐ
end;

// フィルタ係数を取得
procedure GetFilterCoef(Count, SamplingFreq, CutOffFreq:Cardinal; Buffer:PDoubleArray; WindowType:Word);
var
      i,n : Integer;
      Window : TWindowFunc;
      BufPtr : PDouble;
begin
      Window := nil;

      (* フィルタ係数算出 *)

      // メモリ領域チェック
      {if Length(Buffer^) < Count then
      begin
            Exception.Create('Scanty Buffer-Size on GetFilterCoef');
      end;
      }

      BufPtr := PDouble(Buffer);

      // フィルタ係数を書き込む
      // （フーリエ級数打ち切り形）
      for i := 1 to Count do
      begin
            n := i - ((Count-1) div 2);

            if n <> 0 then
                  BufPtr^ := (1 / (pi*n)) * sin(n * (2 * pi * CutOffFreq)/SamplingFreq)
            else
                  BufPtr^ := (2 * CutOffFreq * Count) / SamplingFreq * 0.005;

            BufPtr^ := BufPtr^ * 180;

            Inc(BufPtr);
      end;

      // 窓関数を選ぶ
      case WindowType of
        Win_Kaiser : Window := Window_Kaiser;
        Win_Hanning : Window := Window_Hanning;
        Win_Hamming : Window := Window_Hamming;
        Win_Blackman : Window := Window_Blackman;
        Win_Sin : Window := Window_sin;
        Win_Rectangler : Window := Window_Rectangler;
      else
        Exception.Create('Unknown Window-Type on GetFilterCoef');
      end;

      // 窓関数を適用する
      if Assigned(Window) then Window(Count, Buffer);
end;

// フィルタをかける
function FilteringProcess(Filter:PDoubleArray; FilterSize:Cardinal; InData:TList):Integer;
var
      i : Integer;
      Temp : Extended;
begin
      Temp := 0;

      for i := 0 to FilterSize-1 do
      begin
            if i < InData.Count then
                  Temp := Temp + PDouble(Filter)^ * SmallInt(InData[i]);
                  
            //Temp := Temp + 0.1 * SmallInt(InData[i]);
            Inc(PDouble(Filter));
      end;

      {if InData.Count > 0 then
            Result := Round( Temp / InData.Count )
      else
            Result := 0;
      }

      Result := Round(Temp / Indata.Count);
end;

{procedure LowPassFilterForWaveFile(Src:String; Dest:String; CutOffFreq:Cardinal;
                                          WindowType:Word; ProgressFunc:TProgressFunc = nil);

  procedure AddToListLimit(List:TList; Item:Pointer; LimitCount:Cardinal);
  begin
      if List.Count > LimitCount then List.Delete(0);
      List.Add(Item);
  end;

  function SubPtr(A,B:Pointer):Cardinal; register;
  begin
      Result := Abs(Dword(A) - Dword(B));
  end;

const
      //MaxFilterCount = 205;    // 125 samples
      MaxFilterCount = 175;
      ReadBufferSize = 5242880;
var
      Reader : TFileStream;
      //Reader : TMemoryStream;
      Writer : TFileStream;
      ReadBuffer : PByte;
      ReadBufPtr : PByte;
      ReadBufLimit : PByte;
      ReadBufSz : Cardinal;
      WriteBuffer : TMemoryStream;
      Header : TRawWaveHeader;
      wTmp : SmallInt;
      bTmp : Byte;
      HisL, HisR : TList;
      ChFlag : Boolean;
      ReadSize : Cardinal;
      ReadSum : Cardinal;
      UnitSize : Byte;
      Info : TWaveInfo;
      Filter : Array [1..MaxFilterCount] of Double;
begin
      HisL := TList.Create;
      HisR := TList.Create;
      try
            Info := AnalyzeWaveFile(Src);

            Reader := TFileStream.Create(Src, fmOpenRead);
            //Reader := TMemoryStream.Create;
            Writer := TFileStream.Create(Dest, fmOpenWrite or fmCreate);
            WriteBuffer := TMemoryStream.Create;
            GetMem(ReadBuffer, ReadBufferSize+100);
            try
                  //Reader.LoadFromFile(Src);

                  // set buffer setting
                  ReadBufPtr := ReadBuffer;
                  ReadBufLimit := ReadBuffer; Inc(ReadBufLimit, ReadBufferSize);

                  // init chflag
                  ChFlag := false;

                  Header.RIFFMark := 'RIFF';
                  Header.PreFileSize := Info.FileSize;
                  Header.WaveMark := 'WAVE';
                  Header.FmtMark := 'fmt ';
                  Header.FmtSize := sizeof(Header.FmtChunk);
                  Header.FmtChunk := Info.FormatChuck;
                  Header.DataMark := 'data';
                  Header.WaveSize := Info.DataLength;

                  Reader.Seek(Info.DataPosition, soFromBeginning);
                  ReadSum := Info.DataPosition;

                  UnitSize := (Header.FmtChunk.BitRate div 8);

                  //Reader.

                  // read wave header
                  //Reader.Read(Header, sizeof(Header));

                  // write wave header
                  Writer.Write(Header, sizeof(Header));

                  WriteBuffer.SetSize(2*1024*1024);

                  // get coeff filter
                  GetFilterCoef(
                        MaxFilterCount,
                        Header.FmtChunk.SamplingRate,
                        CutOffFreq,
                        @Filter,
                        WindowType
                  );

                  if Info.IsWaveFile then
                  begin
                        while true do
                        begin
                              if (ReadBufPtr = ReadBufLimit) or (ReadBufPtr = ReadBuffer) then
                              begin
                                    ReadBufSz := Reader.Read(ReadBuffer^, ReadBufferSize);
                                    ReadBufPtr := ReadBuffer;
                                    ReadBufLimit := ReadBuffer; Inc(ReadBufLimit, ReadBufSz);
                                    //ReadBuffer.Clear;
                                    //ReadBuffer.CopyFrom(Reader, 2*1024*1024);
                                    //ReadBuffer.Seek(0,0);
                              end;

                              //Inc(ReadBufPtr);

                              // Bit converter
                              case Header.FmtChunk.BitRate of
                                    8 :
                                    begin
                                          bTmp := ReadBufPtr^;
                                          Inc(ReadBufPtr);

                                          // read 1 sample
                                          //ReadSize := ReadBuffer.Read(bTmp, sizeof(bTmp));
                                          wTmp := (bTmp - 128) * 257;
                                    end;

                                    16:
                                    begin
                                          wTmp := PSmallInt(ReadBufPtr)^;
                                          Inc(PSmallInt(ReadBufPtr));

                                          // read 1 sample
                                          //ReadSize := ReadBuffer.Read(wTmp, sizeof(wTmp));
                                    end;
                              end;

                              if (ReadBufSz <> ReadBufferSize) and
                                    (Dword(ReadBufPtr) > Dword(ReadBufLimit)) then break;

                              //Inc(ReadSum, ReadSize);
                              //if ReadSize = 0 then break;

                              case Header.FmtChunk.Channels of
                                    2 :
                                    begin
                                          if not ChFlag then
                                          begin
                                                // L Ch.
                                                AddToListLimit(HisL, Pointer(wTmp), MaxFilterCount);

                                                // Filtering
                                                wTmp := FilteringProcess(@Filter, MaxFilterCount, HisL);

                                          end else
                                          begin
                                                // R Ch.
                                                AddToListLimit(HisR, Pointer(wTmp), MaxFilterCount);

                                                // Filtering
                                                wTmp := FilteringProcess(@Filter, MaxFilterCount, HisR);
                                          end;

                                          ChFlag := not ChFlag;
                                    end;

                                    1 :
                                    begin
                                          // L Ch.
                                          AddToListLimit(HisL, Pointer(wTmp), MaxFilterCount);

                                          // Filtering
                                          wTmp := FilteringProcess(@Filter, MaxFilterCount, HisL);
                                    end;
                              else
                                    Exception.Create('Unknown Channels on LowPassFilterForWaveFile');
                              end;
                              

                              // Bit converter
                              case Header.FmtChunk.BitRate of
                                    8 :
                                    begin
                                          // Write 1 sample
                                          bTmp := Round(wTmp / 257) + 128;
                                          WriteBuffer.Write(bTmp, sizeof(bTmp));
                                          //Writer.Write(bTmp, sizeof(bTmp));
                                    end;

                                    16:
                                    begin
                                          // Write 1 sample
                                          WriteBuffer.Write(wTmp, sizeof(wTmp));
                                          //Writer.Write(wTmp, sizeof(wTmp));
                                    end;
                              end;

                              // Buffer is 2Mbytes
                              if WriteBuffer.Size > 2097152 then
                              begin
                                    //Writer.Write(PByte(WriteBuffer.Memory)^, WriteBuffer.Position);
                                    WriteBuffer.Seek(0,0);
                                    Writer.CopyFrom(WriteBuffer, WriteBuffer.Size);
                                    WriteBuffer.Clear;
                              end;

                              // Callback
                              if Assigned(ProgressFunc) then
                              begin
                                   ProgressFunc(Reader.Position - ReadBufSz + Subptr(ReadBuffer, ReadBufPtr), Reader.Size);
                              end;

                        end;

                        // 残りバッファを書き込む
                        if WriteBuffer.Size > 0 then
                        begin
                              Writer.Write(PByte(WriteBuffer.Memory)^, WriteBuffer.Position);
                        end;
                  end;

            finally
                  Reader.Free;
                  Writer.Free;
                  WriteBuffer.Free;
                  FreeMem(ReadBuffer);
            end;
      finally
            HisL.Free;
            HisR.Free;
      end;
end;
}

end.