Formant resonator, speech processing 2

feck-lowpass.pl

#!/usr/bin/perl

use warnings;
use strict;
#use GD::Graph;

my $frequency = 0;
my $amplitude = 1;
my $phase = 0;

my $freq_step = 20;
my $f_cutoff = 100;

my (@freqs, @amps, @phases);

my $PI = 3.14159;
my $T_s = 0.0001;
my $F_s = 10000;

my $F_n = 500;
my $B_n = 200;

my $A_co = 0.0609;
my $B_co = 0.9391;
my $B_res = 2 * exp(-$PI * $B_n * $T_s) * cos(2*$PI * $F_n * $T_s);
my $C_res = -exp(-2 * $PI * $B_n * $T_s);
my $A_res = 1 - $B_res - $C_res;


my @x_n = ();
my $x_n_val = 9500;

my @times = ();
my $ctime = 0;

for (my $i = 0; $i <= 5001; $i++) {
	if ($i == 1) {
		$x_n[$i] = $x_n_val;
	} else {
		$x_n[$i] = 0;
	}
	$times[$i] = $ctime;
	$ctime += $T_s;
}

sub log10 {
	my $in = $_[0];
	return log($in) / log(10);
}

my $cur_freq = $frequency;
for (my $i = 0; $i <= 5000; $i++) {
	my $tmp_sq = $cur_freq / $f_cutoff;

	my $amplitude = 1;
	if ($cur_freq != 0) {
		my $amp_sq_a = $A_res / sqrt(2 * $PI * $cur_freq) * (1 - $C_res) - $B_res;
		my $amp_sq_b = sin(2 * $PI * $cur_freq * $T_s) * (1 + $C_res);
		$amplitude = $amp_sq_a * $amp_sq_a + $amp_sq_b * $amp_sq_b;
	}

	my $amp_db = 20*log10($amplitude);

	$phase = 0.0-(atan2($cur_freq, $f_cutoff));
	print "Freq\t$cur_freq\tAmp\t$amplitude\tAmpDB\t$amp_db\tPhase\t$phase\n";

	push @freqs, $cur_freq;
	push @amps, $amplitude;
	push @phases, $phase;
	$cur_freq = $cur_freq + $freq_step;
}

five-formant-feck.pl

#!/usr/bin/perl

use warnings;
use strict;
use utf8;

my $frequency = 0;
my $amplitude = 1;
my $phase = 0;

my $freq_step = 20;
my $f_cutoff = 250;

my (@freqs, @amps, @phases);

sub log10 {
	my $in = $_[0];
	return log($in) / log(10);
}

my $PI = 3.14159;
my $T_s = 0.0001;
my $F_s = 10000;

my $F_n = 500;
my $B_n = 200;


# both F and B are frequencies in Hz.
sub amp_response {
	my $F_n = shift; # resonant frequency
	my $B_n = shift; # bandwidth
	my $f = shift;

	if ($f == 0) {
		return 1;
	}

	my $B_res = 2 * exp(-$PI * $B_n * $T_s) * cos(2 * $PI * $F_n * $T_s);
	my $C_res = -exp(-2 * $PI * $B_n * $T_s);
	my $A_res = 1 - $B_res - $C_res;

	my $amp_sq_a = $A_res / sqrt(2 * $PI * $f) * (1 - $C_res) - $B_res;
	my $amp_sq_b = sin(2 * $PI * $f * $T_s) * (1 + $C_res);
	my $amplitude = $amp_sq_a * $amp_sq_a + $amp_sq_b * $amp_sq_b;
	$amplitude;
}

sub phase_response {
	my $F_n = shift;
	my $B_n = shift;
	my $f = shift;

	my $B_res = 2 * exp(-$PI * $B_n * $T_s) * cos(2 * $PI * $F_n * $T_s);
	my $C_res = -exp(-2 * $PI * $B_n * $T_s);
	my $A_res = 1 - $B_res - $C_res;

	return 0.0 - atan2(sin(2 * $PI * $f * $T_s) * (1 + $C_res), cos(2 * $PI * $f * $T_s) * (1 - $C_res) - $B_res);
}

sub impulse {
	my $impulse = shift;
	my $prev = shift;
	my $pprev = shift;

	my $B_res = 2 * exp(-$PI * $B_n * $T_s) * cos(2 * $PI * $F_n * $T_s);
	my $C_res = -exp(-2 * $PI * $B_n * $T_s);
	my $A_res = 1 - $B_res - $C_res;

	my $imp = $A_res * $impulse + $B_res * $prev + $C_res * $pprev;
	$imp;
}

sub abs_complex {
	my $real = shift;
	my $imaginary = shift;

	return sqrt(($real * $real) + ($imaginary * $imaginary));
}

my $cur_freq = $frequency;
my $prev = 0;
my $pprev = 0;

my $prev2 = 0;
my $pprev2 = 0;
my $prev3 = 0;
my $pprev3 = 0;
my $prev4 = 0;
my $pprev4 = 0;
my $prev5 = 0;
my $pprev5 = 0;

my $cur_imp = 0;

for (my $i = 0; $i <= 5000; $i++) {

	my $F_1 = 500;
	my $B_1 = 100;
	my $F_2 = 1500;
	my $B_2 = 100;
	my $F_3 = 2500;
	my $B_3 = 100;
	my $F_4 = 3500;
	my $B_4 = 100;
	my $F_5 = 4500;
	my $B_5 = 100;
	my $aresp1 = amp_response($F_1, $B_1, $cur_freq);
	my $presp1 = phase_response($F_1, $B_1, $cur_freq);
	my $amp_db1 = 20*log10($aresp1);
	my $aresp2 = amp_response($F_2, $B_2, $cur_freq);
	my $presp2 = phase_response($F_2, $B_2, $cur_freq);
	my $amp_db2 = 20*log10($aresp2);
	my $aresp3 = amp_response($F_3, $B_3, $cur_freq);
	my $presp3 = phase_response($F_3, $B_3, $cur_freq);
	my $amp_db3 = 20*log10($aresp3);
	my $aresp4 = amp_response($F_4, $B_4, $cur_freq);
	my $presp4 = phase_response($F_4, $B_4, $cur_freq);
	my $amp_db4 = 20*log10($aresp4);
	my $aresp5 = amp_response($F_5, $B_5, $cur_freq);
	my $presp5 = phase_response($F_5, $B_5, $cur_freq);
	my $amp_db5 = 20*log10($aresp5);
	my $amp_all = $amp_db1 + $amp_db2 + $amp_db3 + $amp_db4 + $amp_db5;
	my $phase_all = $presp1 + $presp2 + $presp3 + $presp4 + $presp5;
#	print "Amp. Resp 1\t$aresp1\tAmp. Resp 2\t$aresp2\tAmp. Resp 3\t$aresp3\tAmp. Resp 4\t$aresp4\tAmp. Resp 5\t$aresp5\tAmp DB\t$amp_all\tPAll\t$phase_all\n";
	# Need higher pole correction, to correct for missing higher formants (poles)

	#Aco * impulse + Bco  * prev + Cco * pprev
	#Aco * AB4 + Bco  * prev + Cco * pprev

	if ($i == 2) {
		$cur_imp = 10000;
	} else {
		$cur_imp = 0;
	}

	my $imp1 = impulse($cur_imp, $prev, $pprev);
	my $imp2 = impulse($imp1, $prev2, $pprev2);
	my $imp3 = impulse($imp2, $prev3, $pprev3);
	my $imp4 = impulse($imp3, $prev4, $pprev4);
	my $imp5 = impulse($imp4, $prev5, $pprev5);
	$pprev = $prev;
	$prev = $imp1;
	$pprev2 = $prev2;
	$prev2 = $imp2;
	$pprev3 = $prev3;
	$prev3 = $imp3;
	$pprev4 = $prev4;
	$prev4 = $imp4;
	$pprev5 = $prev5;
	$prev5 = $imp5;

	print "Impulse: $imp1 $imp2 $imp3 $imp4 $imp5\n";

	$cur_freq += $freq_step;
}

five-formant-resonator.pl

#!/usr/bin/perl

use warnings;
use strict;
use utf8;

my $frequency = 0;
my $amplitude = 1;
my $phase = 0;

my $freq_step = 10;
my $f_cutoff = 250;

my (@freqs, @amps, @phases);

sub log10 {
	my $in = $_[0];
	return log($in) / log(10);
}

# both F and B are frequencies in Hz.
sub amp_response {
	my $F_n = shift; # resonant frequency
	my $B_n = shift; # bandwidth
	my $f = shift;

	# ∞ ??
	return 0 if ($B_n == 0 && $F_n == 0 && $f == 0);

	my $B_n_sq_div4 = ($B_n * $B_n)/4;
	my $F_n_sq = ($F_n * $F_n);
	my $tmp_to_sq = ($f - $F_n);
	my $tmp_div_a = sqrt($B_n_sq_div4 + ($tmp_to_sq * $tmp_to_sq));
	my $tmp_to_sq2 = ($f + $F_n);
	my $tmp_div_b = sqrt($B_n_sq_div4 + ($tmp_to_sq2 * $tmp_to_sq2));

	my $res = (($F_n_sq + $B_n_sq_div4) / ($tmp_div_a * $tmp_div_b));
	$res;
}

sub phase_response {
	my $F_n = shift;
	my $B_n = shift;
	my $f = shift;

	my $first = atan2((2 * ($f - $F_n)), $B_n);
	my $second = atan2((2 * ($f + $F_n)), $B_n);
	return 0.0 - $first - $second; 
}

sub abs_complex {
	my $real = shift;
	my $imaginary = shift;

	return sqrt(($real * $real) + ($imaginary * $imaginary));
}

my $cur_freq = $frequency;
for (my $i = 0; $i <= 5000; $i++) {

	my $F_1 = 500;
	my $B_1 = 100;
	my $F_2 = 1500;
	my $B_2 = 100;
	my $F_3 = 2500;
	my $B_3 = 100;
	my $F_4 = 3500;
	my $B_4 = 100;
	my $F_5 = 4500;
	my $B_5 = 100;
	my $aresp1 = amp_response($F_1, $B_1, $cur_freq);
	my $presp1 = phase_response($F_1, $B_1, $cur_freq);
	my $amp_db1 = 20*log10($aresp1);
	my $aresp2 = amp_response($F_2, $B_2, $cur_freq);
	my $presp2 = phase_response($F_2, $B_2, $cur_freq);
	my $amp_db2 = 20*log10($aresp2);
	my $aresp3 = amp_response($F_3, $B_3, $cur_freq);
	my $presp3 = phase_response($F_3, $B_3, $cur_freq);
	my $amp_db3 = 20*log10($aresp3);
	my $aresp4 = amp_response($F_4, $B_4, $cur_freq);
	my $presp4 = phase_response($F_4, $B_4, $cur_freq);
	my $amp_db4 = 20*log10($aresp4);
	my $aresp5 = amp_response($F_5, $B_5, $cur_freq);
	my $presp5 = phase_response($F_5, $B_5, $cur_freq);
	my $amp_db5 = 20*log10($aresp5);
	my $amp_all = $amp_db1 + $amp_db2 + $amp_db3 + $amp_db4 + $amp_db5;
	print "Amp. Resp 1\t$aresp1\tAmp. Resp 2\t$aresp2\tAmp. Resp 3\t$aresp3\tAmp. Resp 4\t$aresp4\tAmp. Resp 5\t$aresp5\tAmp DB\t$amp_all\n";
	# Need higher pole correction, to correct for missing higher formants (poles)
	$cur_freq += $freq_step;
}

formant-resonator.pl

#!/usr/bin/perl

use warnings;
use strict;
use utf8;

my $frequency = 0;
my $amplitude = 1;
my $phase = 0;

my $freq_step = 10;
my $f_cutoff = 250;

my (@freqs, @amps, @phases);

sub log10 {
	my $in = $_[0];
	return log($in) / log(10);
}

# both F and B are frequencies in Hz.
sub amp_response {
	my $F_n = shift; # resonant frequency
	my $B_n = shift; # bandwidth
	my $f = shift;

	# ∞ ??
	return 0 if ($B_n == 0 && $F_n == 0 && $f == 0);

	my $B_n_sq_div4 = ($B_n * $B_n)/4;
	my $F_n_sq = ($F_n * $F_n);
	my $tmp_to_sq = ($f - $F_n);
	my $tmp_div_a = sqrt($B_n_sq_div4 + ($tmp_to_sq * $tmp_to_sq));
	my $tmp_to_sq2 = ($f + $F_n);
	my $tmp_div_b = sqrt($B_n_sq_div4 + ($tmp_to_sq2 * $tmp_to_sq2));

	my $res = (($F_n_sq + $B_n_sq_div4) / ($tmp_div_a * $tmp_div_b));
	$res;
}

sub phase_response {
	my $F_n = shift;
	my $B_n = shift;
	my $f = shift;

	my $first = atan2((2 * ($f - $F_n)), $B_n);
	my $second = atan2((2 * ($f + $F_n)), $B_n);
	return 0.0 - $first - $second; 
}

sub abs_complex {
	my $real = shift;
	my $imaginary = shift;

	return sqrt(($real * $real) + ($imaginary * $imaginary));
}

my $cur_freq = $frequency;
for (my $i = 0; $i <= 5000; $i++) {

	my $F_n = 250;
	my $B_n = 50;
	my $aresp = amp_response($F_n, $B_n, $cur_freq);
	my $presp = phase_response($F_n, $B_n, $cur_freq);
	my $amp_db = 20*log10($aresp);
	print "Amp. Resp\t$aresp\tPh. Resp.\t$presp\tAmp DB\t$amp_db\n";
	$cur_freq += $freq_step;
}

z_lowpass.pl

#!/usr/bin/perl

use warnings;
use strict;
use utf8;

my $F_s = 10000;
my $T_s = 0.0001;
my $f_cut = 1;


my (@freqs, @amps, @phases);

sub log10 {
	my $in = $_[0];
	return log($in) / log(10);
}

my @x_n = ();
my @y_n = ();

my $PI = 3.14159;

my $B_coeff = exp(-2*$PI*$f_cut*$T_s);
my $A_coeff = 1 - $B_coeff;

sub do_y_n {
	my $x_n = $_[0];
	my $prev;
	my $prev = $_[1];
	return $A_coeff * $x_n + $B_coeff * $prev;
}

# fill input: x(n) 
for (my $i = 0; $i < 5001; $i++) {
	if ($i == 0) {
		$x_n[$i] = 9500;
	} else {
		$x_n[$i] = 0;
	}
}

print "A coeff\t$A_coeff\tB coeff\t$B_coeff\n\n";

my $last_y = 0;
my $tmp_y = 0;

for (my $i = 0; $i < 5001; $i++) {
	$tmp_y = do_y_n($x_n[$i], $last_y);
	print "y_n\t$tmp_y\n";
	$last_y = $tmp_y;
}