ディジタル信号処理による信号発生とエイリアス(4)

ナイキスト周波数以下の帯域に制限する方法として、実用性は薄いけれど理想的な手段である、項数を制限したフーリエ級数の計算を採用します。
これまでの例では ですから、基本波から 5 倍波までを加算すれば、理想のこぎり波に理想特性のフィルタをかけたことに相当します。
前の「gen_saw1.py」プログラムを実行すると、理想のこぎり波をサンプリングした波形が「saw0.wav」ファイルに、フーリエ級数の計算による波形が「saw1.wav」ファイルに、それぞれ生成されます。
フーリエ級数の計算は、複数の sin 関数の値を求めて加算するわけですから、実用的ではありませんが、理論的にはエイリアス成分をゼロにすることができます。
前のプログラムでは、のこぎり波の１周期がサンプリング周期の整数倍ぴったりの場合でしたが、今回のプログラムでは、のこぎり波の１周期がサンプリング周期の整数倍から、わずかにずれた場合の波形を生成します。
１周期 の波形を 93 周期発生させると 11 × 93 = 1023 サンプルとなります。
ここで、のこぎり波の周期を 11 よりわずかに大きくして、93 周期分が 1024 サンプルとなる周期を選ぶと、スペクトルを求めるのに 1024 点 (以上の) FFT を利用できます。
具体的には、48 [kHz] × 1023 / (11 × 1024) = 4.359375 [kHz] に選びます。
次のプログラムを実行すると、「saw2.wav」ファイルには、左チャンネル：理想のこぎり波のサンプリング波形、右チャンネル：フーリエ級数によるのこぎり波が生成されます。
同時に作成される「saw4.wav」ファイルと、帯域制限の実用的な方法については、次回以降に説明します。

# gen_saw2.py - Sawtooth wave generation
import wave, struct, math;
from math import *;
nperi = 11
fs = 48000 # sampling rate = 48 kHz
maxlevel = 0x6000
nrep = 48
ncyc = 1024 
freq = 1023.0 / ncyc
nharm = (nperi-1) / 2
n_ovs = 2
#
def waveinit():
    """Initialize wave header."""
    ww.setnchannels(2) # stereo
    ww.setsampwidth(2) # 16 bit
    ww.setframerate(fs) 
#    
def saw2(x):
    """Sampled ideal sawtooth wave."""
    return -maxlevel*atan2(cos(2.0*pi*x),
                           sin(2.0*pi*x))/pi
#
def saw3(x):
    """Band limited saw by Fourier series."""
    s = 0.0
    for i in range(1, nharm+1):
        s = s + sin(i * 2.0 * pi * x)/i
    return -maxlevel * 2.0 * s / pi
#
def saw4(x):
    """Nx oversampled/decimated saw2."""
    s = 0.0
    for i in range(n_ovs):
        s = s + saw2(x + i*freq/(nperi*n_ovs))
    return s/n_ovs
#
def waveout(func1, func2):
    """Wave data output."""
    waveinit()
    for i in range(nrep): 
        x = 0.5
        for j in range(ncyc):
            l = func1(x/nperi)
            r = func2(x/nperi)
            lr = struct.pack("hh", l, r)
            ww.writeframesraw(lr)
            x = x + freq
            if (x > nperi) : x = x - nperi;
    ww.close()
# Lch = sampled idal sawtooth wave
# Rch = band limited saw wave by Fourier series
ww = wave.open("saw2.wav", "wb")
waveout(saw2, saw3)
# Lch = Nx oversampled/decimated idal saw wave
# Rch = band limited saw wave by Fourier series
ww = wave.open("saw4.wav", "wb")
waveout(saw4, saw3)

「saw2.wav」ファイルをループ再生して、音声出力をオシロスコープで観察すると、右チャンネルのフーリエ級数による方法の波形は、安定していてジッタがないのに対し、左チャンネルの理想のこぎり波のサンプリングによる方法の波形は、最大１サンプル周期の幅となるジッタが乗っていることが分かります。
音を耳で聞くと、右チャンネルの音は澄んでいますが、左チャンネルには、低い音が混じっていることが分かります。
これは、11 倍波のエイリアスがベースバンド帯域では 46.875 Hz となり、その他のエイリアスと合成されて 46.875 Hz の、のこぎり波となるためです。
そのあたりの事情は、次回以降に説明します。