帰ってきたネコブログ

FIRフィルタとか使う時に必要なので、描いてみた。

use v5.14;
use strict;
use warnings;

use Imager;
use Math::Trig qw/pi/;

use constant STEP_X => 4;
use constant SCALE_Y => 256.0;
use constant N => 101;

my ( $width, $height ) = ( 601, 400 );
my ( $x0, $y0 ) = ( int($width / 2), ($height - 50) );

{
    my $img = Imager->new(
        xsize => $width, ysize => $height, channels => 4 );
    draw_graduation( $img, Imager::Color->new(64,64,64) );

    draw_window( $img, kaiser_window(N, 0.000), 'red' );
    draw_window( $img, kaiser_window(N, 3.395), 'green' );
    draw_window( $img, kaiser_window(N, 5.653), 'orange' );
    draw_window( $img, kaiser_window(N, 7.857), 'yellow' );

    $img->write( file => sprintf('kaiser_%02d.png', N) ) or die $img->errstr;
}

{
    my $img = Imager->new(
        xsize => $width, ysize => $height, channels => 4 );
    draw_graduation( $img, Imager::Color->new(64,64,64) );

    draw_window( $img, hamming_window(N), 'red' );
    draw_window( $img, hann_window(N), 'green' );
    draw_window( $img, kaiser_window(N, 7.857), 'yellow' );

    $img->write( file => sprintf('window_%02d.png', N) ) or die $img->errstr;
}

sub kaiser_window {
    my ( $n, $alpha ) = @_;

    my @ret = map {
        bessel( $alpha * sqrt(1.0 - ($_ ** 2.0)) ) / bessel( $alpha );
    } map {
        ($_ - int($n / 2)) / int($n / 2);
    } 0..($n - 1);

    return \@ret;
}

sub hann_window {
    my $n = shift;

    my @ret = map {
        my $x = 2.0 * pi() * $_;
        0.5 + (0.5 * cos($x));
    } map {
        ($_ - int($n / 2)) / ($n - 1);
    } 0..($n - 1);

    return \@ret;
}

sub hamming_window {
    my $n = shift;

    my @ret = map {
        my $x = 2.0 * pi() * $_;
        0.54 + (0.46 * cos($x));
    } map {
        ($_ - int($n / 2)) / ($n - 1);
    } 0..($n - 1);

    return \@ret;
}

sub bessel {
    my $x = shift;

    my $total = 1;
    my $factorial = 1;
    foreach my $n ( 1..20 ) {
       $factorial *= $n;
       my $tmp = (($x / 2) ** $n) / $factorial;
       $total += ($tmp * $tmp);
    }

    return $total;
}

sub draw_window {
    my ( $img, $data, $color ) = @_;

    my $img_tmp = Imager->new(
        xsize => $img->getwidth(), ysize => $img->getheight(), channels => 4 );

    my $n = scalar( @{$data} );
    for (my $i=0; $i<$n; $i++) {
        my $dx = $i - int($n / 2);
        my $dy = int($data->[$i] * SCALE_Y);

        $img_tmp->line(
            x1 => ($x0 + ($dx * STEP_X)), y1 => $y0,
            x2 => ($x0 + ($dx * STEP_X)), y2 => ($y0 - $dy),
            color => $color );
    }

    $img->compose( src => $img_tmp, opacity => 0.25, combine => 'add' );

    for (my $i=0; $i<$n; $i++) {
        my $dx = $i - int($n / 2);
        my $dy = int($data->[$i] * SCALE_Y);

        my ( $x, $y ) = ( $x0 + ($dx * STEP_X), $y0 - $dy );
        $img->box(
            xmin => $x - 1, ymin => $y - 1, xmax => $x + 1, ymax => $y + 1,
            color => $color, filled => 0 );
    }
}

sub draw_graduation {
    my ( $img, $color ) = @_;

    $img->box( filled => 1, color => 'black' );

    $img->line(
        x1 => $x0, y1 => $y0,
        x2 => $x0, y2 => 0,
        color => $color );
    $img->line(
        x1 => 0, y1 => $y0,
        x2 => ($width - 2), y2 => $y0,
        color => $color );
}

カイザー窓の元ネタはこの本だけど、ハン窓とハミング窓は調べて追加した。

カイザー窓はこんな感じ。

ハミング窓（赤）、ハン窓（緑）、カイザー窓（黄）はこんな感じ。

で、さっきの本だとFIRフィルタの係数を算出できなくて困ったので、
家にある適当な本を使って算出した。
たぶん、sinc関数を使うべきところで、sinc関数に触れてなかったので、
最初はソースコードが載っている入門本が良いと思う。

追記 2016/06/04
ソースコードの59行目と72行目の分母が間違ってたので修正しました。
なので、2枚目のグラフも差し替えました。
— ここまで追記 —

おしまい。

カテゴリー: Perl, コンピューター by ネコ on 2016年04月30日 | タグ: メモ > 音楽