2015/07/30

LADSPA IIR LPF Butterworth 2pole

1poleと同じように2poleのLPFをLADSPAで作ってみる。

上記式から2poleのプロトタイプを計算

LPFを作りたいので、上記式 s に次式を代入する。
LPF: ((z-1)/(z+1)) / (2πfc)

fcの中はこんなかんじ。ωcはカットオフ周波数で、fsはサンプリング周波数。

ブロック図は下図を採用。



LPF 2pole Butterworth ソースコード

処理部分以外は1poleとほとんど同じ。
/* namagi_lpf_2p.c 2015.07.29
 * LPF Butterworth 2pole
 * compile windows
 * gcc -shared -o namagi_lpf_2p.dll namagi_lpf_2p.c -ID
 * 
 * compile Ubuntu
 * gcc -fPIC -DPIC -shared -nostartfiles
   -o namagi_lpf_2p.so namagi_lpf_2p.c
 */
/**********************************************************/
#include <stdlib.h>
#include <string.h>
#include <stdio.h>
#include <math.h>
/**********************************************************/
/*LADSPAのヘッダファイルを読み込む*/
#include "ladspa.h"
/**********************************************************/
/* マイクロ定数を定義 コントロールと出力 */
#define LPF_CONTROL 0
#define LPF_INPUT1  1
#define LPF_OUTPUT1 2
/**********************************************************/
#ifdef WIN32
   int bIsFirstTime = 1;
   void _init();
   #define _WINDOWS_DLL_EXPORT_ __declspec(dllexport)
#else
   #define _WINDOWS_DLL_EXPORT_ 
#endif
/**********************************************************/
/* 構造体で入出力とコントロールを定義 */
typedef struct {
  LADSPA_Data   m_fSampleRate; /* LADSPA_Data は float */
  LADSPA_Data * m_pfControlValue; /* コントロール */
  LADSPA_Data * m_pfInputBuffer1; /* mono 入力 */
  LADSPA_Data * m_pfOutputBuffer1; /* mono 出力 */
} Lpf_t;
/**********************************************************/
/* インスタンス化 */
LADSPA_Handle instantiatelpf(const LADSPA_Descriptor * Descriptor,
       unsigned long SampleRate) {
  Lpf_t * pssLpf;
  pssLpf =  (Lpf_t *)malloc(sizeof(Lpf_t));
  pssLpf -> m_fSampleRate = (LADSPA_Data)SampleRate;
  return; //malloc(sizeof(Lpf_t));
  }
/**********************************************************/
/* portに接続 コントロールと入出力 */
void connectPortTolpf(LADSPA_Handle Instance,
         unsigned long Port,
         LADSPA_Data * DataLocation) 
{
  Lpf_t * pslpf;
  pslpf = (Lpf_t *)Instance;

  switch (Port) {
  case LPF_CONTROL:
    pslpf->m_pfControlValue = DataLocation;
    break;
  case LPF_INPUT1:
    pslpf->m_pfInputBuffer1 = DataLocation;
    break;
  case LPF_OUTPUT1:
    pslpf->m_pfOutputBuffer1 = DataLocation;
    break;
  }
}
/**********************************************************/
/* LPF Butterworth 2pole 処理 */
void runlpf(LADSPA_Handle Instance, unsigned long SampleCount) 
{ 
  LADSPA_Data * pfInput;
  LADSPA_Data * pfOutput;
  LADSPA_Data fwc; /* カットオフ周波数 Hz */
  unsigned long lSampleIndex;
  Lpf_t * pslpf;
  pslpf = (Lpf_t *)Instance;
  
  fwc = *(pslpf    -> m_pfControlValue);
  pfInput = pslpf  -> m_pfInputBuffer1;
  pfOutput = pslpf -> m_pfOutputBuffer1;
  
  float Z1 = 0; /* tap */
  float Z2 = 0;
  float Zn = 0;
  
  /* 配列の作成 */
  float a[3];
  float b[3];
  
  /* 各係数を求める */
  float fc = tan(M_PI * fwc / pslpf -> m_fSampleRate) / (2 * M_PI); 
  float A = 2 * M_PI * fc;
  float P = 1.414;
  
  a[0]=  pow(A, 2)+A*P+1;
  a[1]= -1*(2* pow(A, 2)-2)/a[0];
  a[2]= -1*(pow(A, 2)-A*P+1)/a[0];
  b[0]=  pow(A, 2)/a[0];
  b[1]=  2* pow(A, 2)/a[0];
  b[2]=  pow(A, 2)/a[0];
  a[0]=  1.0;
  
  /* 入出力 */
  for (lSampleIndex = 0; lSampleIndex < SampleCount; lSampleIndex++) {
    *(pfOutput) = ((*(pfInput)+(Z1 * a[1])+(Z2 * a[2])) * b[0])
                  + (Z1 * b[1]) + (Z2 * b[2]);
    Zn = Z1;
    Z1 = *(pfInput) + (Z1 * a[1]) + (Z2 * a[2]);
    Z2 = Zn;
    *(pfInput++);
    *(pfOutput++);
  }
}
/**********************************************************/
/* 開放 */
void cleanuplpf(LADSPA_Handle Instance) {
  free(Instance);
}
/**********************************************************/
/*  */
LADSPA_Descriptor * g_psMonoDescriptor = NULL;
/**********************************************************/
/* 最初に呼び出される部分 インターフェイスや名前などの記述 */
void _init() 
{
  char ** pcPortNames;
  int ports_num = 3;
  LADSPA_PortDescriptor * piPortDescriptors;
  LADSPA_PortRangeHint  * psPortRangeHints;
  
    g_psMonoDescriptor
      = (LADSPA_Descriptor *)malloc(sizeof(LADSPA_Descriptor));
    g_psMonoDescriptor -> UniqueID
      = 1048;
    g_psMonoDescriptor -> Label
      = strdup("LPF_2pole");
    g_psMonoDescriptor -> Properties
      = LADSPA_PROPERTY_HARD_RT_CAPABLE;
    g_psMonoDescriptor -> Name 
      = strdup("Namagi: LPF 2Pole ver.150729");
    g_psMonoDescriptor -> Maker
      = strdup("Namagi Products");
    g_psMonoDescriptor -> Copyright
      = strdup("None");
    g_psMonoDescriptor -> PortCount
      = ports_num;
    piPortDescriptors
      = (LADSPA_PortDescriptor *)calloc(ports_num, 
        sizeof(LADSPA_PortDescriptor));
    g_psMonoDescriptor->PortDescriptors
      = (const LADSPA_PortDescriptor *)piPortDescriptors;
    piPortDescriptors[LPF_CONTROL]
      = LADSPA_PORT_INPUT | LADSPA_PORT_CONTROL;
    piPortDescriptors[LPF_INPUT1]
      = LADSPA_PORT_INPUT | LADSPA_PORT_AUDIO;
    piPortDescriptors[LPF_OUTPUT1]
      = LADSPA_PORT_OUTPUT | LADSPA_PORT_AUDIO;
    pcPortNames
      = (char **)calloc(ports_num, sizeof(char *));
    g_psMonoDescriptor->PortNames 
      = (const char **)pcPortNames;
    pcPortNames[LPF_CONTROL]
      = strdup("wc(Hz)");
    pcPortNames[LPF_INPUT1]
      = strdup("Input");
    pcPortNames[LPF_OUTPUT1]
      = strdup("Output");
    psPortRangeHints
      = ((LADSPA_PortRangeHint *)
        calloc(ports_num, sizeof(LADSPA_PortRangeHint)));
    g_psMonoDescriptor->PortRangeHints
      = (const LADSPA_PortRangeHint *)psPortRangeHints;
    psPortRangeHints[LPF_CONTROL].HintDescriptor
      = (LADSPA_HINT_BOUNDED_BELOW 
        | LADSPA_HINT_LOGARITHMIC | LADSPA_HINT_DEFAULT_1);

    psPortRangeHints[LPF_CONTROL].HintDescriptor
      = (LADSPA_HINT_BOUNDED_BELOW | 
      LADSPA_HINT_BOUNDED_ABOVE | LADSPA_HINT_DEFAULT_LOW);
    psPortRangeHints[LPF_CONTROL].LowerBound 
      = 0;
    psPortRangeHints[LPF_CONTROL].UpperBound
      = 20000;
    psPortRangeHints[LPF_INPUT1].HintDescriptor
      = 0;
    psPortRangeHints[LPF_OUTPUT1].HintDescriptor
      = 0;
    g_psMonoDescriptor -> instantiate 
      = instantiatelpf;
    g_psMonoDescriptor -> connect_port 
      = connectPortTolpf;
    g_psMonoDescriptor -> activate
      = NULL;
    g_psMonoDescriptor -> run
      = runlpf;
    g_psMonoDescriptor -> run_adding
      = NULL;
    g_psMonoDescriptor -> set_run_adding_gain
      = NULL;
    g_psMonoDescriptor -> deactivate
      = NULL;
    g_psMonoDescriptor -> cleanup
      = cleanuplpf;
}
/**********************************************************/
/* 開放 インターフェイスまわり */
void deleteDescriptor(LADSPA_Descriptor * psDescriptor) {
  unsigned long lIndex;
  if (psDescriptor) {
    free((char *)psDescriptor->Label);
    free((char *)psDescriptor->Name);
    free((char *)psDescriptor->Maker);
    free((char *)psDescriptor->Copyright);
    free((LADSPA_PortDescriptor *)psDescriptor->PortDescriptors);
    for (lIndex = 0; lIndex < psDescriptor->PortCount; lIndex++)
      free((char *)(psDescriptor->PortNames[lIndex]));
    free((char **)psDescriptor->PortNames);
    free((LADSPA_PortRangeHint *)psDescriptor->PortRangeHints);
    free(psDescriptor);
  }
}
/**********************************************************/
void _fini() {
  deleteDescriptor(g_psMonoDescriptor);
}
/**********************************************************/
_WINDOWS_DLL_EXPORT_
const LADSPA_Descriptor * ladspa_descriptor(unsigned long Index) {
#ifdef WIN32
  if (bIsFirstTime) {
  _init(); 
  bIsFirstTime = 0;
  }
#endif
  if (Index == 0)
    return g_psMonoDescriptor;
  else
    return NULL;
}

Audacity(Win)で実行すると以下のような画面になる。

上段がホワイトノイズにで、下段がLPF 2pole カットオフ周波数1000Hzを適用したところ。

周波数スペクトル。1poleに比べて急なカーブを描く。Audacity標準のLow Pass Filterとほとんど同じ結果を出すので、間違ってはいなさそう。



sound programming 目次はこちら