mkunfilter
This program has been superceded by an equivalent shell script in the mkfilter distribution by the same name.
mkunfilter is a simple Perl script to create an impulse response that has the inverse frequency response of the one given. It was originally created for a technique in room impulse response measurement.
Usage:
perl mkunfilter.pl original.wav reverse.wav 500 4
The third argument is the output length in samples, and the fourth is the maximum gain on any given frequency band. Due to noisy inputs, precision issues, and the fact that frequencies can drop out, it is vitally important that the maximum gain be set to a reasonable value.
#!/usr/bin/perl
#############################################################################
# mkunfilter.pl - Create a reverse FIR filter for the given impulse response
# The latest version of this can be found at:
# http://encryptio.com/code/mkunfilter
#############################################################################
# Copyright (c) 2009 Chris Kastorff
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are
# met:
#
# * Redistributions of source code must retain the above copyright
# notice, this list of conditions, and the following disclaimer.
# * Redistributions in binary form must reproduce the above copyright
# notice, this list of conditions and the following disclaimer in the
# documentation and/or other materials provided with the distribution.
# * The name Chris Kastorff may not be used to endorse or promote
# products derived from this software without specific prior written
# permission.
#############################################################################
use warnings;
use strict;
our $VERSION = 0.2;
use Math::FFT;
use Audio::SndFile;
use Math::Trig qw/ pi /;
use POSIX qw/ ceil floor /;
use List::Util qw/ max /;
my ($inputpath, $outputpath, $outputlength, $forcelimit) = @ARGV;
die unless $inputpath and $outputpath and $outputlength and $forcelimit;
my $lf = Audio::SndFile->open("<", $inputpath) or die;
die if $lf->channels != 1;
my $sf = Audio::SndFile->open(">", $outputpath,
type => "wav", # TODO: autodetect
endianness => $lf->endianness,
samplerate => $lf->samplerate,
subtype => "pcm_24",
channels => $lf->channels
);
my @data = $lf->unpackf_double($lf->frames);
my $wantsize = 2 ** ceil(log(scalar @data)/log(2));
push @data, 0 while @data < $wantsize;
# take fft
@data = map +($_,0), @data; # to complex numbers
my $infft = Math::FFT->new(\@data)->cdft;
# extract frequency response curve
my @responsecurve = ();
for my $i ( 0 .. floor($#$infft/2) ) {
push @responsecurve, sqrt $infft->[$i*2]**2 + $infft->[$i*2+1]**2;
}
# invert frequency response curve
for my $i ( 0 .. $#responsecurve ) {
if ( $responsecurve[$i] == 0 ) {
$responsecurve[$i] = $forcelimit;
} else {
$responsecurve[$i] = 1/$responsecurve[$i];
$responsecurve[$i] = $forcelimit if $responsecurve[$i] > $forcelimit;
}
}
@responsecurve = map $_/$forcelimit, @responsecurve;
# take inverse fft
my $coeff = [map +($_,0), @responsecurve];
my $outfft = Math::FFT->new($coeff)->invcdft($coeff);
# build final impulse response
@data = map $outfft->[$_*2], 0 .. floor($#$outfft/2); # extract real part
unshift @data, pop @data for 1 .. $outputlength/2; # shift
@data = @data[0..$outputlength-1]; # truncate
# hamming window
for my $i ( 0 .. $#data ) {
$data[$i] *= 0.54 - 0.46 * cos(2*pi*$i/@data);
}
# write out
$sf->packf_double(@data);
__END__
Revision history:
0.2:
* Complex arithmetic optimization/simplification, 10x speed increase
0.1:
* First Release
#############################################################################
# mkunfilter.pl - Create a reverse FIR filter for the given impulse response
# The latest version of this can be found at:
# http://encryptio.com/code/mkunfilter
#############################################################################
# Copyright (c) 2009 Chris Kastorff
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are
# met:
#
# * Redistributions of source code must retain the above copyright
# notice, this list of conditions, and the following disclaimer.
# * Redistributions in binary form must reproduce the above copyright
# notice, this list of conditions and the following disclaimer in the
# documentation and/or other materials provided with the distribution.
# * The name Chris Kastorff may not be used to endorse or promote
# products derived from this software without specific prior written
# permission.
#############################################################################
use warnings;
use strict;
our $VERSION = 0.2;
use Math::FFT;
use Audio::SndFile;
use Math::Trig qw/ pi /;
use POSIX qw/ ceil floor /;
use List::Util qw/ max /;
my ($inputpath, $outputpath, $outputlength, $forcelimit) = @ARGV;
die unless $inputpath and $outputpath and $outputlength and $forcelimit;
my $lf = Audio::SndFile->open("<", $inputpath) or die;
die if $lf->channels != 1;
my $sf = Audio::SndFile->open(">", $outputpath,
type => "wav", # TODO: autodetect
endianness => $lf->endianness,
samplerate => $lf->samplerate,
subtype => "pcm_24",
channels => $lf->channels
);
my @data = $lf->unpackf_double($lf->frames);
my $wantsize = 2 ** ceil(log(scalar @data)/log(2));
push @data, 0 while @data < $wantsize;
# take fft
@data = map +($_,0), @data; # to complex numbers
my $infft = Math::FFT->new(\@data)->cdft;
# extract frequency response curve
my @responsecurve = ();
for my $i ( 0 .. floor($#$infft/2) ) {
push @responsecurve, sqrt $infft->[$i*2]**2 + $infft->[$i*2+1]**2;
}
# invert frequency response curve
for my $i ( 0 .. $#responsecurve ) {
if ( $responsecurve[$i] == 0 ) {
$responsecurve[$i] = $forcelimit;
} else {
$responsecurve[$i] = 1/$responsecurve[$i];
$responsecurve[$i] = $forcelimit if $responsecurve[$i] > $forcelimit;
}
}
@responsecurve = map $_/$forcelimit, @responsecurve;
# take inverse fft
my $coeff = [map +($_,0), @responsecurve];
my $outfft = Math::FFT->new($coeff)->invcdft($coeff);
# build final impulse response
@data = map $outfft->[$_*2], 0 .. floor($#$outfft/2); # extract real part
unshift @data, pop @data for 1 .. $outputlength/2; # shift
@data = @data[0..$outputlength-1]; # truncate
# hamming window
for my $i ( 0 .. $#data ) {
$data[$i] *= 0.54 - 0.46 * cos(2*pi*$i/@data);
}
# write out
$sf->packf_double(@data);
__END__
Revision history:
0.2:
* Complex arithmetic optimization/simplification, 10x speed increase
0.1:
* First Release
