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

前回の、位相をズラした理想のこぎり波をサンプリングした波形を左右のチャンネルに納めた wave ファイルを生成する Python プログラムを作りました。
後日説明する方法による、位相ズレがちゃんと現れるような手段で作成した波形の wave ファイルも同時に生成されます。
処理系としては、Windows 上の Python 2.5.1 を使いました。
言語として「Python」を採用したのは、wave ファイルを簡単に扱えるモジュールが標準で付属しているためで、Python であることに特別な意味はありません。
今回、初めて Python を使ったので、Python らしい書き方ができているわけでもありません。

# gen_saw1.py - Sawtooth wave generation
import wave, struct, math
from math import *
nperi = 11
ofs = 0.25
fs = 48000 # sampling rate = 48 kHz
maxlevel = 0x6000
ncyc = fs / nperi # cycles per 1 sec
nharm = (nperi-1) / 2
def waveinit():
    """Initialize wave header."""
    ww.setnchannels(2) # stereo
    ww.setsampwidth(2) # 16 bit
    ww.setframerate(fs)
#    
def saw0(x):
    """Sampled ideal sawtooth wave."""
    return -maxlevel*atan2(cos(2.0*pi*x),
                           sin(2.0*pi*x))/pi
#
def saw1(x):
    """Band limited saw wave."""
    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 waveout(func1, func2, ofs):
    """Wave data output."""
    waveinit()
    for j in range(ncyc):
        for i in range(nperi):
            x = float(i) + 0.5
            l = func1(x/nperi)
            r = func2((x+ofs)/nperi)
            lr = struct.pack("hh", l, r)
            ww.writeframesraw(lr)                         
    ww.close()
# sampled idal sawtooth wave
ww = wave.open("saw0.wav", "wb")
waveout(saw0, saw0, ofs)
# band limited saw wave by Fourier series
ww = wave.open("saw1.wav", "wb")
waveout(saw1, saw1, ofs)

L チャンネルと R チャンネルの間に 1/4 サンプル分の位相差を設けてあります。
このプログラムを実行すると、理想のこぎり波をサンプルした場合の「saw0.wav」と、位相差が現れる方法による「saw1.wav」ができます。
両方とも、再生時間が約 1 秒のファイルで、耳で聞くのが目的ではありません。
ループ再生が可能な wave ファイル・プレイヤーで再生して、音声出力波形をオシロスコープで観測するか、あるいは、波形編集ソフトでファイルの中身の波形を見るためのものです。
下の図が波形編集ソフトでみた「saw0.wav」、つまり、理想のこぎり波をサンプリングした波形で、位相差が現れていないことが分かります。

赤の横線が水平軸ですが、水平軸上に乗っているサンプルを比べてみると、わずかにレベルが違うことが分かります。
下の図は、位相差が表れる方法の「saw1.wav」の波形です。

黄色の縦の破線で示されているカーソル位置の補間波形を見ると、L/R 間で位相が違うのが、はっきり分かります。