頼りないニモニック

はっきりいって個人の日記レベル

オーバーサンプリングによる折り返し雑音の除去

(この記事は以前に http://nanatomo.com/program/cs/707 で公開されていたものです。)

音の折り返し雑音の除去に関する日本語ページが無いため、このページを記述しました。

結論から言えば、オーバーサンプリングによる折り返し雑音(エイリアス)除去は可能です。ただし処理がどうしても重くなってしまうため、除去しながら再生するというのは難しいと思われます。

このページでは折り返し雑音の原理やローパスフィルタの仕組みなどは記述しません。

除去手順

折り返し雑音が含まれている原信号から雑音を除去する手順は以下の通りです。

  1. アップサンプリング ・・・ 原信号のサンプリング周波数fsをn倍に変換する (nは整数が望ましい)
  2. フィルタ処理・・ カットオフ周波数 fc を fc < fs/2 として決定し、ローパスフィルタにかける。
  3. ダウンサンプリング ・・・ 原信号の fs にサンプリング周波数を戻す。

以上の処理を行うことによって除去が可能です。サンプリング周波数変換に用いる n は整数のほうが単純な処理で済みます。またフィルタ処理における fc はフィルタの実装により具体的な値は変わります。

テスト環境ではフィルタ処理にFIR型のローパスフィルタを適用しました。FIR型はIIR型に比べると処理コストはかかりますが、よりよいフィルタリング結果を得られます。

結果

n を増やしていくことによって折り返し雑音が除去されていく様子を周波数スペクトルでまとめました。図中の fp はフィルタ処理時に用いたサンプリング周波数で、 fp = fsn です。

使用波形: ノコギリ波 (500Hz)

f:id:nanase_t:20130412162116p:plain

まとめ

n を増やしていくことによって折り返し雑音は限りなく除去することが可能です。ただしそれに従って全体の処理やメモリ消費が比例して増加するため、n を不用意に増加させるのは得策ではありません。

また n をあまり増加させなくとも、(原信号にもよるが)人間の耳では認識されない程度に折り返し雑音は除去できます。

ソースコード

ライブラリとして ALSharp が必要です。

using System;
using ALSharp;

namespace Antialiasing
{
  class Program
  {
    static void Main(string[] args)
    {
      OpenAL.AlutInit();

      const int Factor = 16;
      const int BufferSize = 2048;
      const double SamplingFreq = 44100;
      const double PI2 = Math.PI * 2.0;
      const double PI2_SR = Math.PI * 2.0 / (SamplingFreq * Factor);
      const double Freq = 500;

      double phase = 0;

      FirFilter f = new FirFilter(SamplingFreq * Factor, 1500, BufferSize * Factor);
      f.SetFilter(FilterType.LowPass, 21000, 0);

      double[] input = new double[BufferSize * Factor];
      double[] output = new double[BufferSize * Factor];
      Random r = new Random();

      Func<double[], int, int, int> process = (b, o, c) =>
      {
        for (int i = 0; i < BufferSize * Factor; i++)
          input[i] = (((PI2_SR * Freq * (i + phase) + Math.PI) % PI2) / Math.PI - 1.0) * 0.8;

        f.Filtering(input, output);

        for (int i = 0; i < BufferSize; i++)
          b[i] = output[i * Factor];

        phase += BufferSize * Factor;

        return c;
      };

      using (RawWavePlayer p = new RawWavePlayer(process, new PlayerSettings() { ChannelCount = 1, SamplingFrequency = (int)SamplingFreq, BufferCount = 16, BufferSize = BufferSize * 2 }))
        Console.ReadLine();

      OpenAL.AlutExit();
    }
  }

  class FirFilter
  {
    private readonly double samplingFreq;
    private readonly double delta;
    private readonly double[] filter, x;
    private readonly int delayer, size;

    public FirFilter(double samplingFreq, double delta, int bufferSize)
    {
      delta /= samplingFreq;

      this.samplingFreq = samplingFreq;
      this.delta = delta;
      this.size = bufferSize;

      this.delayer = (int)(3.1 / delta + 0.5) - 1;
      if (delayer % 2 == 1)
        delayer++;

      filter = new double[delayer + 1];
      x = new double[this.size + this.delayer];
    }

    public void SetFilter(FilterType type, double fe1, double fe2)
    {
      int offset = this.delayer / 2;

      fe1 /= this.samplingFreq;
      fe2 /= this.samplingFreq;

      double fe1_2 = 2.0 * fe1,
           fe2_2 = 2.0 * fe2,
           fe1_PI2 = 2.0 * Math.PI * fe1,
           fe2_PI2 = 2.0 * Math.PI * fe2;

      switch (type)
      {
        case FilterType.LowPass:
          for (int i = -this.delayer / 2, length = this.delayer / 2; i <= length; i++)
            this.filter[offset + i] = fe1_2 * Sinc(fe1_PI2 * i);

          break;
        case FilterType.HighPass:
          for (int i = -this.delayer / 2, length = this.delayer / 2; i <= length; i++)
            this.filter[offset + i] = Sinc(Math.PI * i) - fe1_2 * Sinc(fe1_PI2 * i);

          break;
        case FilterType.BandPass:
          for (int i = -this.delayer / 2, length = this.delayer / 2; i <= length; i++)
            this.filter[offset + i] = fe2_2 * Sinc(fe2_PI2 * i)
                        - fe1_2 * Sinc(fe1_PI2 * i);

          break;
        case FilterType.BandElimination:
          for (int i = -this.delayer / 2, length = this.delayer / 2; i <= length; i++)
            this.filter[offset + i] = Sinc(Math.PI * i)
                        - fe2_2 * Sinc(fe2_PI2 * i)
                        + fe1_2 * Sinc(fe1_PI2 * i);

          break;
        default:
          break;
      }

      FirFilter.GetHanningWindow(filter);
    }

    public void Filtering(double[] input, double[] output)
    {
      for (int i = this.delayer, length = this.size + this.delayer; i < length; i++)
        this.x[i] = input[i - this.delayer];

      for (int i = 0; i < this.size; i++)
      {
        output[i] = 0.0;
        for (int j = 0; j <= this.delayer; j++)
          output[i] += this.filter[j] * this.x[this.delayer + i - j];
      }

      for (int i = this.size - this.delayer, j = 0; i < this.size; i++, j++)
        this.x[j] = input[i];
    }

    private static double Sinc(double x)
    {
      return x == 0.0 ? 1.0 : Math.Sin(x) / x;
    }

    private static void GetHanningWindow(double[] array)
    {
      double k = 2.0 * Math.PI / (double)array.Length;

      if (array.Length % 2 == 0)
        for (int i = 0, length = array.Length; i < length; i++)
          array[i] *= 0.5 - 0.5 * Math.Cos(k * i);
      else
        for (int i = 0, length = array.Length; i < length; i++)
          array[i] *= 0.5 - 0.5 * Math.Cos(k * (i + 0.5));
    }
  }

  enum FilterType
  {
    LowPass,
    HighPass,
    BandPass,
    BandElimination
  }
}