kenmirrors
/
freeswitch
mirror of https://github.com/signalwire/freeswitch.git
1
0
Fork 0
Code Issues Projects Releases Wiki Activity

git status -u, learn something new every day.

This commit is contained in:
Brian West 2012-12-20 20:17:20 -06:00
parent 59205c7678
commit 68913681a4
106 changed files with 10920 additions and 0 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

199
libs/libcodec2/README_fdmdv.txt Normal file
View File

@ -0,0 +1,199 @@
README_fdmdv.txt
David Rowe
Created March 2012
Introduction
------------
A 1400 bit/s Frequency Division Multiplexed Digital Voice (FDMDV) modem
based on [1]. Used for digital audio over HF SSB.
The FDMDV modem was first implemented in GNU Octave, then ported to C.
Algorithm development is generally easier in Octave, but for real time
work we need the C version. Automated units tests ensure the
operation of the Octave and C versions are identical.
Quickstart
----------
$ cd codec2-dev
$ ./configure && make
$ cd src
1. Generate some test bits and modulate them:
$ ./fdmdv_get_test_bits test.c2 1400
$ ./fdmdv_mod test.c2 test.raw
$ play -r 8000 -s -2 test.raw
2. Two seconds of test frame data modulated and sent out of sound device:
$ ./fdmdv_get_test_bits - 2800 | ./fdmdv_mod - - | play -t raw -r 8000 -s -2 -
3. Send 14000 modulated bits (10 seconds) to the demod and count errors:
$ ./fdmdv_get_test_bits - 14000 | ./fdmdv_mod - - | ./fdmdv_demod - - demod_dump.txt | ./fdmdv_put_test_bits -
Use Octave to look at plots of 1 second (1400 bits) of modem operation:
$ cd ../octave
$ octave
octave:1> fdmdv_demod_c("../src/demod_dump.txt",1400)
4. Run Octave simulation of entire modem and AWGN channel:
$ cd ../octave
$ octave
octave:1> fdmdv_ut
5. NOTE: If you would like to play modem samples over the air please
convert the 8 kHz samples to 48 kHz. Many PC sound cards have
wildly inaccurate sample clock rates when set to 8 kHz, but seem to
perform OK when set for 48 kHz. If playing and recording files you
can use the sox utility:
$ sox -r 8000 -s -2 modem_sample_8kHz.raw -r 48000 modem_sample_48kHz.wav
For real time applications, the fdmdv.[ch] library includes functions to
convert between 48 and 8 kHz sample rates.
References
----------
[1] http://n1su.com/fdmdv/FDMDV_Docs_Rel_1_4b.pdf
[2] http://n1su.com/fdmdv/
[3] http://www.rowetel.com/blog/?p=2433 "Testing a FDMDV Modem"
[4] http://www.rowetel.com/blog/?p=2458 "FDMDV Modem Page" on David's web site
C Code
------
src/fdmdv_mod.c - C version of modulator that takes a file of bits and
converts it to a raw file of modulated samples.
src/fdmdv_demod.c - C version of demodulator that takes a raw file of
modulated samples and outputs a file of bits.
Optionally dumps demod states to a text file which
can be plotted using the Octave script
fdmdv_demod_c.m
src/fdmdv.h - Header file that exposes FDMDV C API functions. Include
this file in your application program.
src/fdmdv.c - C functions that implement the FDMDV modem.
src/fdmdv-internal.h - internal states and constants for FDMDV modem,
shouldn't be exposed to application program.
unittest/tfdmdv.c - Used to conjunction with unittest/tfdmdv.m to
automatically test C FDMDV functions against
Octave versions.
Octave Scripts
--------------
(Note these require some Octave packages to be installed, see
octave/README.txt).
fdmdv.m - Functions and variables that implement the Octave version of
the FDMDV modem.
fdmdv_ut.m - Unit test for fdmdv Octave code, useful while
developing algorithm. Includes tx/rx plus basic channel
simulation.
Typical run:
octave:6> fdmdv_ut
Eb/No (meas): 7.30 (8.29) dB
bits........: 2464
errors......: 20
BER.........: 0.0081
PAPR........: 13.54 dB
SNR.........: 4.0 dB
It also outputs lots of nice plots that show the
operation of the modem.
For a 1400 bit/s DQPSK modem we expect about 1% BER for
Eb/No = 7.3dB, which corresponds to SNR = 4dB (3kHz
noise BW). The extra dB of measured power is due to the
DBPSK pilot. Currently the noise generation code
doesn't take the pilot power into account, so in this
example the real SNR is actually 5dB.
fdmdv_mod.m - Octave version of modulator that outputs a raw file.
The modulator is driven by a test frame of bits. This
can then be played over a real channel or through a
channel simulator like PathSim. The sample rate can be
changed using "sox" to simulate differences in tx/rx
sample clocks.
To generate 10 seconds of modulated signal:
octave:8> fdmdv_mod("test.raw",1400*10);
fdmdv_demod.m - Demodulator program that takes a raw file as input,
and works out the bit error rate using the known test
frame. Can be used to test the demod performs with
off-air signals, or signals that have been passed
through a channel simulator.
To demodulate 2 seconds of the test.raw file generated
above:
octave:9> fdmdv_demod("test.raw",1400*2);
2464 bits 0 errors BER: 0.0000
It also produces several plots showing the internal
states of the demod. Useful for debugging and
observing what happens with various channels.
fdmdv_demod_c.m - Takes an output text file from the C demod
fdmdv_demod.c and produces plots and measures BER.
Useful for evaluating fdmdv_demod.c performance.
The plots produced are identical to the Octave
version fdmdv_demod.m, allowing direct comparison of
the C and Octave versions.
tfdmdv.m - Automatic tests that compare the Octave and C versions of
the FDMDV modem functions. First run unittest/tfdmdv, this
will generate a text file with test vectors from the C
version. Then run the Octave script tfdmdv and it will
generate Octave versions of the test vectors and compare
each vector with the C equivalent. Its plots the vectors
and and errors (green). Its also produces an automatic
check list based on test results. If the Octave or C modem
code is changed, this script should be used to ensure the
C and Octave versions remain identical.
Modelling sample clock errors using sox
---------------------------------------
This introduces a simulated 1000ppm error:
sox -r 8000 -s -2 mod_dqpsk.raw -s -2 mod_dqpsk_8008hz.raw rate -h 8008
TODO
----
[ ] implement ppm measurements in fdmdv_get_demod_stats()
[ ] try interfering sine wave
+ maybe swept
+ does modem fall over?
[ ] try non-flat channel, e.g. 3dB difference between hi and low tones
+ make sure all estimators keep working
[ ] test rx level sensitivity, i.e. 0 to 20dB attenuation
[ ] make fine freq indep of amplitude
+ use angle rather than imag coord
[ ] document use of fdmdv_ut and fdmdv_demod + PathSim
[ ] more positive form of sync reqd for DV frames?
+ like using coarse_fine==1 to decode valid DV frame bit?
+ when should we start decoding?
[ ] more robust track/acquite state machine?
+ e.g. hang on thru the fades?
[ ] PAPR idea
+ automatically tweak phases to reduce PAPR, e.g. slow variations in freq...
[ ] why is pilot noise_est twice as big as other carriers

109
libs/libcodec2/asterisk/README Normal file
View File

@ -0,0 +1,109 @@
README for codec2/asterisk
Asterisk Codec 2 support
Test Configuration
------------------
Codec 2 is used to trunk calls between two Asterisk boxes:
A - SIP phone - Asterisk A - Codec2 - Asterisk B - SIP Phone - B
The two SIP phones are configured for mulaw.
Building
---------
Asterisk must be patched so that the core understand Codec 2 frames.
1/ First install Codec 2:
david@cool:~$ svn co https://freetel.svn.sourceforge.net/svnroot/freetel/codec2-dev codec2-dev
david@cool:~/codec2-dev$ cd codec2-dev
david@cool:~/codec2-dev$ ./configure && make && sudo make install
david@bear:~/codec2-dev$ sudo ldconfig -v
david@cool:~/codec2-dev$ cd ~
2/ Then build Asterisk with Codec 2 support:
david@cool:~$ tar xvzf asterisk-1.8.9.0.tar.gz
david@cool:~/asterisk-1.8.9.0$ patch -p4 < ~/codec2-dev/asterisk/asterisk-codec2.patch
david@cool:~/asterisk-1.8.9.0$ cp ~/codec2-dev/asterisk/codec_codec2.c .
david@cool:~/asterisk-1.8.9.0$ cp ~/codec2-dev/asterisk/ex_codec2.h ./codecs
david@cool:~/asterisk-1.8.9.0$ ./configure && make ASTLDFLAGS=-lcodec2
david@cool:~/asterisk-1.8.9.0$ sudo make install
david@cool:~/asterisk-1.8.9.0$ sudo make samples
3/ Add this to the end of sip.conf on Asterisk A:
[6013]
type=friend
context=default
host=dynamic
user=6013
secret=6013
canreinvite=no
callerid=6013
disallow=all
allow=ulaw
[potato]
type=peer
username=potato
fromuser=potato
secret=password
context=default
disallow=all
dtmfmode=rfc2833
callerid=server
canreinvite=no
host=cool
allow=codec2
3/ Add this to the end of sip.conf on Asterisk B:
[6014]
type=friend
context=default
host=dynamic
user=6014
secret=6014
canreinvite=no
callerid=6014
disallow=all
allow=ulaw
[potato]
type=peer
username=potato
fromuser=potato
secret=password
context=default
disallow=all
dtmfmode=rfc2833
callerid=server
canreinvite=no
host=bear
allow=codec2
4/ Here is the [default] section of extensions.conf on Asterisk B:
[default]
exten => 6013,1,Dial(SIP/potato/6013)
;
; By default we include the demo. In a production system, you
; probably don't want to have the demo there.
;
;include => demo
5/ After booting see if the codec2_codec2.so module is loaded with "core show translate"
6/ To make a test call dial 6013 on the SIP phone connected to Asterisk B
7/ If codec_codec2.so won't load and you see "can't find codec2_create" try:
david@cool:~/asterisk-1.8.9.0$ touch codecs/codec_codec2.c
david@cool:~/asterisk-1.8.9.0$ make ASTLDFLAGS=-lcodec2
david@cool:~/asterisk-1.8.9.0$ sudo cp codecs/codec_codec2.so /usr/lib/asterisk/modules
david@cool:~/asterisk-1.8.9.0$ sudo asterisk -vvvcn

68
libs/libcodec2/asterisk/asterisk-codec2.patch Normal file
View File

@ -0,0 +1,68 @@
--- /home/david/asterisk-1.8.9.0-orig/include/asterisk/frame.h 2011-12-23 05:08:46.000000000 +1030
+++ /home/david/asterisk-1.8.9.0-codec2/include/asterisk/frame.h 2012-03-27 13:16:55.623452431 +1030
@@ -299,6 +299,7 @@
#define AST_FORMAT_G719 (1ULL << 32)
/*! SpeeX Wideband (16kHz) Free Compression */
#define AST_FORMAT_SPEEX16 (1ULL << 33)
+#define AST_FORMAT_CODEC2 (1ULL << 34)
/*! Raw mu-law data (G.711) */
#define AST_FORMAT_TESTLAW (1ULL << 47)
/*! Reserved bit - do not use */
--- /home/david/asterisk-1.8.9.0-orig/main/frame.c 2010-06-18 02:53:43.000000000 +0930
+++ /home/david/asterisk-1.8.9.0-codec2/main/frame.c 2012-03-28 15:16:16.975581316 +1030
@@ -102,6 +102,7 @@
{ AST_FORMAT_ADPCM, "adpcm" , 8000, "ADPCM", 40, 10, 300, 10, 20 }, /*!< codec_adpcm.c */
{ AST_FORMAT_SLINEAR, "slin", 8000, "16 bit Signed Linear PCM", 160, 10, 70, 10, 20, AST_SMOOTHER_FLAG_BE }, /*!< Signed linear */
{ AST_FORMAT_LPC10, "lpc10", 8000, "LPC10", 7, 20, 20, 20, 20 }, /*!< codec_lpc10.c */
+ { AST_FORMAT_CODEC2, "codec2", 8000, "Codec 2", 7, 20, 20, 20, 20 }, /*!< codec_codec2.c */
{ AST_FORMAT_G729A, "g729", 8000, "G.729A", 10, 10, 230, 10, 20, AST_SMOOTHER_FLAG_G729 }, /*!< Binary commercial distribution */
{ AST_FORMAT_SPEEX, "speex", 8000, "SpeeX", 10, 10, 60, 10, 20 }, /*!< codec_speex.c */
{ AST_FORMAT_SPEEX16, "speex16", 16000, "SpeeX 16khz", 10, 10, 60, 10, 20 }, /*!< codec_speex.c */
@@ -1475,6 +1476,9 @@
samples = 22 * 8;
samples += (((char *)(f->data.ptr))[7] & 0x1) * 8;
break;
+ case AST_FORMAT_CODEC2:
+ samples = 160 * (f->datalen / 7);
+ break;
case AST_FORMAT_ULAW:
case AST_FORMAT_ALAW:
case AST_FORMAT_TESTLAW:
@@ -1519,6 +1523,9 @@
case AST_FORMAT_GSM:
len = (samples / 160) * 33;
break;
+ case AST_FORMAT_CODEC2:
+ len = (samples / 160) * 7;
+ break;
case AST_FORMAT_G729A:
len = samples / 8;
break;
--- /home/david/asterisk-1.8.9.0-orig/main/channel.c 2011-12-17 10:21:13.000000000 +1030
+++ /home/david/asterisk-1.8.9.0-codec2/main/channel.c 2012-03-28 15:05:22.395293391 +1030
@@ -1075,6 +1075,7 @@
/*! Ick, LPC10 sounds terrible, but at least we have code for it, if you're tacky enough
to use it */
AST_FORMAT_LPC10,
+ AST_FORMAT_CODEC2,
/*! G.729a is faster than 723 and slightly less expensive */
AST_FORMAT_G729A,
/*! Down to G.723.1 which is proprietary but at least designed for voice */
--- /home/david/asterisk-1.8.9.0-orig/main/rtp_engine.c 2011-12-30 01:43:03.000000000 +1030
+++ /home/david/asterisk-1.8.9.0-codec2/main/rtp_engine.c 2012-03-28 16:42:02.880872891 +1030
@@ -101,6 +101,7 @@
{{1, AST_FORMAT_SLINEAR}, "audio", "L16", 8000},
{{1, AST_FORMAT_SLINEAR16}, "audio", "L16", 16000},
{{1, AST_FORMAT_LPC10}, "audio", "LPC", 8000},
+ {{1, AST_FORMAT_CODEC2}, "audio", "CODEC2", 8000},
{{1, AST_FORMAT_G729A}, "audio", "G729", 8000},
{{1, AST_FORMAT_G729A}, "audio", "G729A", 8000},
{{1, AST_FORMAT_G729A}, "audio", "G.729", 8000},
@@ -178,6 +179,7 @@
[117] = {1, AST_FORMAT_SPEEX16},
[118] = {1, AST_FORMAT_SLINEAR16}, /* 16 Khz signed linear */
[121] = {0, AST_RTP_CISCO_DTMF}, /* Must be type 121 */
+ [121] = {1, AST_FORMAT_CODEC2},
};
int ast_rtp_engine_register2(struct ast_rtp_engine *engine, struct ast_module *module)

185
libs/libcodec2/asterisk/codec_codec2.c Normal file
View File

@ -0,0 +1,185 @@
/*
* Codec 2 module for Asterisk.
*
* Credit: codec_gsm.c used as a starting point.
*
* Copyright (C) 2012 Ed W and David Rowe
*
* This program is free software, distributed under the terms of
* the GNU General Public License Version 2. See the LICENSE file
* at the top of the source tree.
*/
/*! \file
*
* \brief Translate between signed linear and Codec 2
*
* \ingroup codecs
*/
/*** MODULEINFO
<support_level>core</support_level>
***/
#include "asterisk.h"
#include "asterisk/translate.h"
#include "asterisk/config.h"
#include "asterisk/module.h"
#include "asterisk/utils.h"
#include <codec2.h>
#define BUFFER_SAMPLES 8000
#define CODEC2_SAMPLES 160
#define CODEC2_FRAME_LEN 7
/* Sample frame data */
#include "asterisk/slin.h"
#include "ex_codec2.h"
struct codec2_translator_pvt { /* both codec2tolin and codec2togsm */
struct CODEC2 *codec2;
short buf[BUFFER_SAMPLES]; /* lintocodec2, temporary storage */
};
static int codec2_new(struct ast_trans_pvt *pvt)
{
struct codec2_translator_pvt *tmp = pvt->pvt;
tmp->codec2 = codec2_create(CODEC2_MODE_2400);
return 0;
}
/*! \brief decode and store in outbuf. */
static int codec2tolin_framein(struct ast_trans_pvt *pvt, struct ast_frame *f)
{
struct codec2_translator_pvt *tmp = pvt->pvt;
int x;
int16_t *dst = pvt->outbuf.i16;
int flen = CODEC2_FRAME_LEN;
for (x=0; x < f->datalen; x += flen) {
unsigned char *src;
int len;
len = CODEC2_SAMPLES;
src = f->data.ptr + x;
codec2_decode(tmp->codec2, dst + pvt->samples, src);
pvt->samples += CODEC2_SAMPLES;
pvt->datalen += 2 * CODEC2_SAMPLES;
}
return 0;
}
/*! \brief store samples into working buffer for later decode */
static int lintocodec2_framein(struct ast_trans_pvt *pvt, struct ast_frame *f)
{
struct codec2_translator_pvt *tmp = pvt->pvt;
if (pvt->samples + f->samples > BUFFER_SAMPLES) {
ast_log(LOG_WARNING, "Out of buffer space\n");
return -1;
}
memcpy(tmp->buf + pvt->samples, f->data.ptr, f->datalen);
pvt->samples += f->samples;
return 0;
}
/*! \brief encode and produce a frame */
static struct ast_frame *lintocodec2_frameout(struct ast_trans_pvt *pvt)
{
struct codec2_translator_pvt *tmp = pvt->pvt;
int datalen = 0;
int samples = 0;
/* We can't work on anything less than a frame in size */
if (pvt->samples < CODEC2_SAMPLES)
return NULL;
while (pvt->samples >= CODEC2_SAMPLES) {
/* Encode a frame of data */
codec2_encode(tmp->codec2, (unsigned char*)(pvt->outbuf.c + datalen), tmp->buf + samples);
datalen += CODEC2_FRAME_LEN;
samples += CODEC2_SAMPLES;
pvt->samples -= CODEC2_SAMPLES;
}
/* Move the data at the end of the buffer to the front */
if (pvt->samples)
memmove(tmp->buf, tmp->buf + samples, pvt->samples * 2);
return ast_trans_frameout(pvt, datalen, samples);
}
static void codec2_destroy_stuff(struct ast_trans_pvt *pvt)
{
struct codec2_translator_pvt *tmp = pvt->pvt;
if (tmp->codec2)
codec2_destroy(tmp->codec2);
}
static struct ast_translator codec2tolin = {
.name = "codec2tolin",
.srcfmt = AST_FORMAT_CODEC2,
.dstfmt = AST_FORMAT_SLINEAR,
.newpvt = codec2_new,
.framein = codec2tolin_framein,
.destroy = codec2_destroy_stuff,
.sample = codec2_sample,
.buffer_samples = BUFFER_SAMPLES,
.buf_size = BUFFER_SAMPLES * 2,
.desc_size = sizeof (struct codec2_translator_pvt ),
};
static struct ast_translator lintocodec2 = {
.name = "lintocodec2",
.srcfmt = AST_FORMAT_SLINEAR,
.dstfmt = AST_FORMAT_CODEC2,
.newpvt = codec2_new,
.framein = lintocodec2_framein,
.frameout = lintocodec2_frameout,
.destroy = codec2_destroy_stuff,
.sample = slin8_sample,
.desc_size = sizeof (struct codec2_translator_pvt ),
.buf_size = (BUFFER_SAMPLES * CODEC2_FRAME_LEN + CODEC2_SAMPLES - 1)/CODEC2_SAMPLES,
};
/*! \brief standard module glue */
static int reload(void)
{
return AST_MODULE_LOAD_SUCCESS;
}
static int unload_module(void)
{
int res;
res = ast_unregister_translator(&lintocodec2);
if (!res)
res = ast_unregister_translator(&codec2tolin);
return res;
}
static int load_module(void)
{
int res;
res = ast_register_translator(&codec2tolin);
if (!res)
res=ast_register_translator(&lintocodec2);
else
ast_unregister_translator(&codec2tolin);
if (res)
return AST_MODULE_LOAD_FAILURE;
return AST_MODULE_LOAD_SUCCESS;
}
AST_MODULE_INFO(ASTERISK_GPL_KEY, AST_MODFLAG_DEFAULT, "Codec 2 Coder/Decoder",
.load = load_module,
.unload = unload_module,
.reload = reload,
);

28
libs/libcodec2/asterisk/ex_codec2.h Normal file
View File

@ -0,0 +1,28 @@
/*! \file
* \brief 8-bit raw data
*
* Copyright (C) 2012, 2012 Ed W and David Rowe
*
* Distributed under the terms of the GNU General Public License
*
*/
static uint8_t ex_codec2[] = {
0x3e,0x06,0x4a,0xbb,0x9e,0x40,0xc0
};
static struct ast_frame *codec2_sample(void)
{
static struct ast_frame f = {
.frametype = AST_FRAME_VOICE,
.subclass.codec = AST_FORMAT_CODEC2,
.datalen = sizeof(ex_codec2),
.samples = CODEC2_SAMPLES,
.mallocd = 0,
.offset = 0,
.src = __PRETTY_FUNCTION__,
.data.ptr = ex_codec2,
};
return &f;
}

11
libs/libcodec2/asterisk/make_asterisk_patch.sh Executable file
View File

@ -0,0 +1,11 @@
#!/bin/sh
# Create patch for Codec 2 support inside Asterisk
ORIG=~/asterisk-1.8.9.0-orig
CODEC2=~/asterisk-1.8.9.0-codec2
diff -ruN $ORIG/include/asterisk/frame.h $CODEC2/include/asterisk/frame.h > asterisk-codec2.patch
diff -ruN $ORIG/main/frame.c $CODEC2/main/frame.c >> asterisk-codec2.patch
diff -ruN $ORIG/main/channel.c $CODEC2/main/channel.c >> asterisk-codec2.patch
diff -ruN $ORIG/main/rtp_engine.c $CODEC2/main/rtp_engine.c >> asterisk-codec2.patch

19
libs/libcodec2/fltk/Makefile Normal file
View File

@ -0,0 +1,19 @@
# Requires FLTK 1.3 & Portaudio V19
FLTK_VER = $(shell fltk-config --api-version)
ifneq ($(FLTK_VER),1.3)
$(error Must use FLTK version 1.3, you have $(FLTK_VER))
endif
FLTK_CFLAGS += $(shell fltk-config --ldstaticflags)
CFLAGS = -O3 -g -Wall
LIBS = ../src/.libs/libcodec2.a -lm -lrt -lportaudio -pthread
LC2POC_C = fl_fdmdv.cxx
all: fl_fdmdv
fl_fdmdv: Makefile $(LC2POC_C)
g++ $(LC2POC_C) -I../src/ -o fl_fdmdv $(CFLAGS) $(FLTK_CFLAGS) $(LIBS)
clean:
rm -f fl_fdmdv

1158
libs/libcodec2/fltk/fl_fdmdv.cxx Normal file
View File

File diff suppressed because it is too large Load Diff

17
libs/libcodec2/octave/README.txt Normal file
View File

@ -0,0 +1,17 @@
README.txt
For codec2/octave directory
Created 24 June 2012 by David Rowe
1/ To support some of the Octave scripts (in particular fdmdv) in this
directory the following Octave packages need to be installed:
control image miscellaneous optim signal specfun struct
2/ Download these tar balls from:
http://octave.sourceforge.net/packages.php
3/ Install each package from the Octave command line with:
octave:3> pkg install package_file_name.tar.gz

43
libs/libcodec2/octave/av_imp.m Normal file
View File

@ -0,0 +1,43 @@
% av_imp.m
% David Rowe Aug 2012
% Averages the impulse response samples
function imp = av_imp(imp_filename, period_in_secs, st, en)
f = fopen(imp_filename,"rb");
s = fread(f, Inf, "short")';
Fs = 8000;
n = period_in_secs * Fs;
[r c] = size(s);
imp = zeros(1,n);
for i=1:n:c-n
imp = imp + s(i:i+n-1);
endfor
% user supplies start and end samples after viweing plot
if (nargin == 4)
imp = imp(st:en);
end
% normalise
imp /= sqrt(sum(imp .^ 2));
[h w] = freqz(imp, 1, 4000);
figure(1);
clf;
plot(imp);
figure(2);
clf;
subplot(211)
plot(10*log10(abs(h)))
subplot(212)
plot(angle(h))
endfunction

98
libs/libcodec2/octave/cbphase.m Normal file
View File

@ -0,0 +1,98 @@
% cbphase.m
% David Rowe Aug 2012
% Used to experiment with critical band phase perception and smoothing
function cbphase
Wo = 100.0*pi/4000;
L = floor(pi/Wo);
A = zeros(1,L);
phi = zeros(1,L);
% three harmonics in this band
b = 4; a = b-1; c = b+1;
% set up phases and mags for 2nd order system (see phasesecord.m)
wres = b*Wo;
phi(a) = 3*pi/4 + wres;
phi(b) = pi/2 + wres;
phi(c) = pi/4 + wres;
A(a) = 0.707;
A(b) = 1;
A(c) = 0.707;
% add linear component
phi(1) = pi;
phi(2:L) = phi(2:L) + (2:L)*phi(1);
phi = phi - 2*pi*(floor(phi/(2*pi)) + 0.5);
N = 16000;
Nplot = 250;
s = zeros(1,N);
for m=a:c
s_m = A(m)*cos(m*Wo*(0:(N-1)) + phi(m));
s = s + s_m;
endfor
figure(2);
clf;
subplot(211)
plot((1:L)*Wo*4000/pi, A,'+');
subplot(212)
plot((1:L)*Wo*4000/pi, phi,'+');
%v = A(a)*exp(j*phi(a)) + A(b)*exp(j*phi(b)) + A(c)*exp(j*phi(c));
%compass(v,"r")
%hold off;
% est phi1
diff = phi(b) - phi(a)
sumi = sin(diff);
sumr = cos(diff);
diff = phi(c) - phi(b)
sumi += sin(diff);
sumr += cos(diff);
phi1_ = atan2(sumi, sumr)
s_v = cos(Wo*(0:(N-1)) + phi1_);
figure(1);
clf;
subplot(211)
plot(s(1:Nplot));
hold on;
plot(s_v(1:Nplot),"r");
hold off;
% build (hopefully) perceptually similar phase
phi_(a) = a*phi1_;
phi_(b) = b*phi1_;
phi_(c) = c*phi1_;
s_ = zeros(1,N);
for m=a:c
s_m = A(m)*cos(m*Wo*(0:(N-1)) + phi_(m));
s_ = s_ + s_m;
endfor
subplot(212)
plot(s_(1:Nplot));
gain = 8000;
fs=fopen("orig_ph.raw","wb");
fwrite(fs,gain*s,"short");
fclose(fs);
fs=fopen("mod_ph.raw","wb");
fwrite(fs,gain*s_,"short");
fclose(fs);
endfunction

108
libs/libcodec2/octave/codec2_demo.m Normal file
View File

@ -0,0 +1,108 @@
% Copyright David Rowe 2012
% This program is distributed under the terms of the GNU General Public License
% Version 2
%
% codec2_demo.m
% Designed as an educational tool to explain the operation of Codec 2
% for conference and user group presentations on a projector. An
% alternative to static overhead slides.
%
% Derived from codec2-dev/octave/plamp.m
%
% usage:
% octave:1> plamp("../src/hts2a",40)
%
% Then press:
% c - to cycle through the wavform being displayed on the figure
% n - next frame
% b - back one frame
%
% tip: hold down n or b to animate the display
%
% The text files used as input are generated using c2sim:
%
% /codec2-dev/src$ c2sim ../raw/hts2a.raw --dump hts2a
%
% The Codec 2 README explains how to build c2sim with dump files
% enabled.
function codec2_demo(samname, f)
sn_name = strcat(samname,"_sn.txt");
Sn = load(sn_name);
sw_name = strcat(samname,"_sw.txt");
Sw = load(sw_name);
model_name = strcat(samname,"_model.txt");
model = load(model_name);
figure(1);
k = ' ';
wf = "Sn";
do
if strcmp(wf,"Sn")
clf;
s = [ Sn(2*f-1,:) Sn(2*f,:) ];
plot(s);
axis([1 length(s) -20000 20000]);
end
if (strcmp(wf,"Sw"))
clf;
plot((0:255)*4000/256, Sw(f,:),";Sw;");
end
if strcmp(wf,"SwAm")
Wo = model(f,1);
L = model(f,2);
Am = model(f,3:(L+2));
plot((0:255)*4000/256, Sw(f,:),";Sw;");
hold on;
plot((1:L)*Wo*4000/pi, 20*log10(Am),"+;Am;r");
axis([1 4000 -10 80]);
hold off;
end
if strcmp(wf,"Am")
Wo = model(f,1);
L = model(f,2);
Am = model(f,3:(L+2));
plot((1:L)*Wo*4000/pi, 20*log10(Am),"+;Am;r");
axis([1 4000 -10 80]);
end
% interactive menu
printf("\rframe: %d menu: n-next b-back w-cycle window q-quit", f);
fflush(stdout);
k = kbhit();
if (k == 'n')
f = f + 1;
end
if (k == 'b')
f = f - 1;
end
if (k == 'w')
if strcmp(wf,"Sn")
next_wf = "Sw";
end
if strcmp(wf,"Sw")
next_wf = "SwAm";
end
if strcmp(wf,"SwAm")
next_wf = "Am";
end
if strcmp(wf,"Am")
next_wf = "Sn";
end
wf = next_wf;
end
until (k == 'q')
printf("\n");
endfunction

54
libs/libcodec2/octave/cspec.m Normal file
View File

@ -0,0 +1,54 @@
% cspec.m
% David Rowe Aug 2012
% Used to compare spectromgrams while experimenting with phase
function cspec(s1,s2)
f1 = fopen(s1,"rb");
s1 = fread(f1,Inf,"short");
f2 = fopen(s2,"rb");
s2 = fread(f2,Inf,"short");
Fs = 8000;
spec_win = 512;
state = 's1';
do
if strcmp(state,'s1')
spec(s1,Fs,spec_win);
%title(s1);
end
if strcmp(state,'s2')
spec(s2,Fs,spec_win);
%title(s2);
end
if strcmp(state,'diff')
spec(s1-s2,Fs,spec_win);
%title("difference");
end
printf("\rstate: %s space-toggle d-diff q-quit", state);
fflush(stdout);
k = kbhit();
if k == ' '
if strcmp(state,"diff")
next_state = 's1';
end
if strcmp(state,"s1")
next_state = 's2';
end
if strcmp(state,'s2')
next_state = 's1';
end
end
if k == 'd'
next_state = 'diff';
end
state = next_state;
until (k == 'q')
printf("\n");
endfunction

935
libs/libcodec2/octave/fdmdv.m Normal file
View File

@ -0,0 +1,935 @@
% fdmdv.m
%
% Functions that implement a Frequency Divison Multiplexed Modem for
% Digital Voice (FDMDV) over HF channels.
%
% Copyright David Rowe 2012
% This program is distributed under the terms of the GNU General Public License
% Version 2
%
% reqd to mak sure we get same random bits at mod and demod
rand('state',1);
randn('state',1);
% Constants
global Fs = 8000; % sample rate in Hz
global T = 1/Fs; % sample period in seconds
global Rs = 50; % symbol rate in Hz
global Nc = 14; % number of carriers
global Nb = 2; % Bits/symbol for QPSK modulation
global Rb = Nc*Rs*Nb; % bit rate
global M = Fs/Rs; % oversampling factor
global Nsym = 6; % number of symbols to filter over
global Fsep = 75; % Separation between carriers (Hz)
global Fcentre = 1200; % Centre frequency, Nc/2 carriers below this, N/c carriers above (Hz)
global Nt = 5; % number of symbols we estimate timing over
global P = 4; % oversample factor used for rx symbol filtering
global Nfilter = Nsym*M;
global Nfiltertiming = M+Nfilter+M;
alpha = 0.5;
global snr_coeff = 0.9 % SNR est averaging filter coeff
% root raised cosine (Root Nyquist) filter
global gt_alpha5_root;
gt_alpha5_root = gen_rn_coeffs(alpha, T, Rs, Nsym, M);
% Functions ----------------------------------------------------
% generate Nc+1 QPSK symbols from vector of (1,Nc*Nb) input bits. The
% Nc+1 symbol is the +1 -1 +1 .... BPSK sync carrier
function tx_symbols = bits_to_qpsk(prev_tx_symbols, tx_bits, modulation)
global Nc;
global Nb;
global pilot_bit;
% re-arrange as (Nc,Nb) matrix
tx_bits_matrix = zeros(Nc,Nb);
tx_bits_matrix(1:Nc,1) = tx_bits(1:Nb:Nb*Nc);
tx_bits_matrix(1:Nc,2) = tx_bits(2:Nb:Nb*Nc);
if (strcmp(modulation,'dqpsk'))
% map to (Nc,1) DQPSK symbols
for c=1:Nc
msb = tx_bits_matrix(c,1); lsb = tx_bits_matrix(c,2);
if ((msb == 0) && (lsb == 0))
tx_symbols(c) = prev_tx_symbols(c);
endif
if ((msb == 0) && (lsb == 1))
tx_symbols(c) = j*prev_tx_symbols(c);
endif
if ((msb == 1) && (lsb == 0))
tx_symbols(c) = -prev_tx_symbols(c);
endif
if ((msb == 1) && (lsb == 1))
tx_symbols(c) = -j*prev_tx_symbols(c);
endif
end
else
% QPSK mapping
tx_symbols = -1 + 2*tx_bits_matrix(:,1) - j + 2j*tx_bits_matrix(:,2);
endif
% +1 -1 +1 -1 BPSK sync carrier, once filtered becomes two spectral
% lines at +/- Rs/2
if pilot_bit
tx_symbols(Nc+1) = -prev_tx_symbols(Nc+1);
else
tx_symbols(Nc+1) = prev_tx_symbols(Nc+1);
end
if pilot_bit
pilot_bit = 0;
else
pilot_bit = 1;
end
endfunction
% Given Nc*Nb bits construct M samples (1 symbol) of Nc filtered
% symbols streams
function tx_baseband = tx_filter(tx_symbols)
global Nc;
global M;
global tx_filter_memory;
global Nfilter;
global gt_alpha5_root;
tx_baseband = zeros(Nc+1,M);
% tx filter each symbol, generate M filtered output samples for each symbol.
% Efficient polyphase filter techniques used as tx_filter_memory is sparse
tx_filter_memory(:,Nfilter) = sqrt(2)/2*tx_symbols;
for i=1:M
tx_baseband(:,i) = M*tx_filter_memory(:,M:M:Nfilter) * gt_alpha5_root(M-i+1:M:Nfilter)';
end
tx_filter_memory(:,1:Nfilter-M) = tx_filter_memory(:,M+1:Nfilter);
tx_filter_memory(:,Nfilter-M+1:Nfilter) = zeros(Nc+1,M);
endfunction
% Construct FDM signal by frequency shifting each filtered symbol
% stream. Returns complex signal so we can apply frequency offsets
% easily.
function tx_fdm = fdm_upconvert(tx_filt)
global Fs;
global M;
global Nc;
global Fsep;
global phase_tx;
global freq;
tx_fdm = zeros(1,M);
% Nc/2 tones below centre freq
for c=1:Nc/2
for i=1:M
phase_tx(c) = phase_tx(c) * freq(c);
tx_fdm(i) = tx_fdm(i) + tx_filt(c,i)*phase_tx(c);
end
end
% Nc/2 tones above centre freq
for c=Nc/2+1:Nc
for i=1:M
phase_tx(c) = phase_tx(c) * freq(c);
tx_fdm(i) = tx_fdm(i) + tx_filt(c,i)*phase_tx(c);
end
end
% add centre pilot tone
c = Nc+1;
for i=1:M
phase_tx(c) = phase_tx(c) * freq(c);
pilot(i) = 2*tx_filt(c,i)*phase_tx(c);
tx_fdm(i) = tx_fdm(i) + pilot(i);
end
% Scale such that total Carrier power C of real(tx_fdm) = Nc. This
% excludes the power of the pilot tone.
% We return the complex (single sided) signal to make frequency
% shifting for the purpose of testing easier
tx_fdm = 2*tx_fdm;
endfunction
% Frequency shift each modem carrier down to Nc+1 baseband signals
function rx_baseband = fdm_downconvert(rx_fdm, nin)
global Fs;
global M;
global Nc;
global Fsep;
global phase_rx;
global freq;
rx_baseband = zeros(1,nin);
% Nc/2 tones below centre freq
for c=1:Nc/2
for i=1:nin
phase_rx(c) = phase_rx(c) * freq(c);
rx_baseband(c,i) = rx_fdm(i)*phase_rx(c)';
end
end
% Nc/2 tones above centre freq
for c=Nc/2+1:Nc
for i=1:nin
phase_rx(c) = phase_rx(c) * freq(c);
rx_baseband(c,i) = rx_fdm(i)*phase_rx(c)';
end
end
% Pilot
c = Nc+1;
for i=1:nin
phase_rx(c) = phase_rx(c) * freq(c);
rx_baseband(c,i) = rx_fdm(i)*phase_rx(c)';
end
endfunction
% Receive filter each baseband signal at oversample rate P
function rx_filt = rx_filter(rx_baseband, nin)
global Nc;
global M;
global P;
global rx_filter_memory;
global Nfilter;
global gt_alpha5_root;
global Fsep;
rx_filt = zeros(Nc+1,nin*P/M);
% rx filter each symbol, generate P filtered output samples for each symbol.
% Note we keep memory at rate M, it's just the filter output at rate P
N=M/P;
j=1;
for i=1:N:nin
rx_filter_memory(:,Nfilter-N+1:Nfilter) = rx_baseband(:,i:i-1+N);
rx_filt(:,j) = rx_filter_memory * gt_alpha5_root';
rx_filter_memory(:,1:Nfilter-N) = rx_filter_memory(:,1+N:Nfilter);
j+=1;
end
endfunction
% LPF and peak pick part of freq est, put in a function as we call it twice
function [foff imax pilot_lpf S] = lpf_peak_pick(pilot_baseband, pilot_lpf, nin)
global M;
global Npilotlpf;
global Npilotcoeff;
global Fs;
global Mpilotfft;
global pilot_coeff;
% LPF cutoff 200Hz, so we can handle max +/- 200 Hz freq offset
pilot_lpf(1:Npilotlpf-nin) = pilot_lpf(nin+1:Npilotlpf);
j = 1;
for i = Npilotlpf-nin+1:Npilotlpf
pilot_lpf(i) = pilot_baseband(j:j+Npilotcoeff-1) * pilot_coeff';
j++;
end
% decimate to improve DFT resolution, window and DFT
Mpilot = Fs/(2*200); % calc decimation rate given new sample rate is twice LPF freq
h = hanning(Npilotlpf);
s = pilot_lpf(1:Mpilot:Npilotlpf) .* h(1:Mpilot:Npilotlpf)';
s = [s zeros(1,Mpilotfft-Npilotlpf/Mpilot)];
S = fft(s, Mpilotfft);
% peak pick and convert to Hz
[imax ix] = max(S);
r = 2*200/Mpilotfft; % maps FFT bin to frequency in Hz
if ix > Mpilotfft/2
foff = (ix - Mpilotfft - 1)*r;
else
foff = (ix - 1)*r;
endif
endfunction
% Estimate frequency offset of FDM signal using BPSK pilot. This is quite
% sensitive to pilot tone level wrt other carriers
function [foff S1 S2] = rx_est_freq_offset(rx_fdm, pilot, pilot_prev, nin)
global M;
global Npilotbaseband;
global pilot_baseband1;
global pilot_baseband2;
global pilot_lpf1;
global pilot_lpf2;
% down convert latest nin samples of pilot by multiplying by
% ideal BPSK pilot signal we have generated locally. This
% peak of the resulting signal is sensitive to the time shift between
% the received and local version of the pilot, so we do it twice at
% different time shifts and choose the maximum.
pilot_baseband1(1:Npilotbaseband-nin) = pilot_baseband1(nin+1:Npilotbaseband);
pilot_baseband2(1:Npilotbaseband-nin) = pilot_baseband2(nin+1:Npilotbaseband);
for i=1:nin
pilot_baseband1(Npilotbaseband-nin+i) = rx_fdm(i) * conj(pilot(i));
pilot_baseband2(Npilotbaseband-nin+i) = rx_fdm(i) * conj(pilot_prev(i));
end
[foff1 max1 pilot_lpf1 S1] = lpf_peak_pick(pilot_baseband1, pilot_lpf1, nin);
[foff2 max2 pilot_lpf2 S2] = lpf_peak_pick(pilot_baseband2, pilot_lpf2, nin);
if max1 > max2
foff = foff1;
else
foff = foff2;
end
endfunction
% Estimate optimum timing offset, re-filter receive symbols
function [rx_symbols rx_timing env] = rx_est_timing(rx_filt, rx_baseband, nin)
global M;
global Nt;
global Nc;
global rx_filter_mem_timing;
global rx_baseband_mem_timing;
global P;
global Nfilter;
global Nfiltertiming;
global gt_alpha5_root;
% nin adjust
% --------------------------------
% 120 -1 (one less rate P sample)
% 160 0 (nominal)
% 200 1 (one more rate P sample)
adjust = P - nin*P/M;
% update buffer of Nt rate P filtered symbols
rx_filter_mem_timing(:,1:(Nt-1)*P+adjust) = rx_filter_mem_timing(:,P+1-adjust:Nt*P);
rx_filter_mem_timing(:,(Nt-1)*P+1+adjust:Nt*P) = rx_filt(:,:);
% sum envelopes of all carriers
env = sum(abs(rx_filter_mem_timing(:,:))); % use all Nc+1 carriers for timing
%env = abs(rx_filter_mem_timing(Nc+1,:)); % just use BPSK pilot
[n m] = size(env);
% The envelope has a frequency component at the symbol rate. The
% phase of this frequency component indicates the timing. So work out
% single DFT at frequency 2*pi/P
x = env * exp(-j*2*pi*(0:m-1)/P)';
% map phase to estimated optimum timing instant at rate M
% the M/4 part was adjusted by experiment, I know not why....
rx_timing = angle(x)*M/(2*pi) + M/4;
if (rx_timing > M)
rx_timing -= M;
end
if (rx_timing < -M)
rx_timing += M;
end
% rx_baseband_mem_timing contains M + Nfilter + M samples of the
% baseband signal at rate M this enables us to resample the filtered
% rx symbol with M sample precision once we have rx_timing
rx_baseband_mem_timing(:,1:Nfiltertiming-nin) = rx_baseband_mem_timing(:,nin+1:Nfiltertiming);
rx_baseband_mem_timing(:,Nfiltertiming-nin+1:Nfiltertiming) = rx_baseband;
% sample right in the middle of the timing estimator window, by filtering
% at rate M
s = round(rx_timing) + M;
rx_symbols = rx_baseband_mem_timing(:,s+1:s+Nfilter) * gt_alpha5_root';
endfunction
% Phase estimation function - probably won't work over a HF channel
% Tries to operate over a single symbol but uses phase information from
% all Nc carriers which should increase the SNR of phase estimate.
% Maybe phase is coherent over a couple of symbols in HF channel,not
% sure but it's worth 3dB so worth experimenting or using coherent as
% an option.
function rx_phase = rx_est_phase(prev_rx_symbols, rx_symbols)
% modulation strip
rx_phase = angle(sum(rx_symbols .^ 4))/4;
endfunction
% convert symbols back to an array of bits
function [rx_bits sync_bit f_err phase_difference] = qpsk_to_bits(prev_rx_symbols, rx_symbols, modulation)
global Nc;
global Nb;
global Nb;
if (strcmp(modulation,'dqpsk'))
% extra 45 degree clockwise lets us use real and imag axis as
% decision boundaries
phase_difference = zeros(Nc+1,1);
phase_difference(1:Nc) = rx_symbols(1:Nc) .* conj(prev_rx_symbols(1:Nc)) * exp(j*pi/4);
% map (Nc,1) DQPSK symbols back into an (1,Nc*Nb) array of bits
for c=1:Nc
d = phase_difference(c);
if ((real(d) >= 0) && (imag(d) >= 0))
msb = 0; lsb = 0;
endif
if ((real(d) < 0) && (imag(d) >= 0))
msb = 0; lsb = 1;
endif
if ((real(d) < 0) && (imag(d) < 0))
msb = 1; lsb = 0;
endif
if ((real(d) >= 0) && (imag(d) < 0))
msb = 1; lsb = 1;
endif
rx_bits(2*(c-1)+1) = msb;
rx_bits(2*(c-1)+2) = lsb;
end
% Extract DBPSK encoded Sync bit
phase_difference(Nc+1,1) = rx_symbols(Nc+1) .* conj(prev_rx_symbols(Nc+1));
if (real(phase_difference(Nc+1)) < 0)
sync_bit = 1;
f_err = imag(phase_difference(Nc+1));
else
sync_bit = 0;
f_err = -imag(phase_difference(Nc+1));
end
% pilot carrier gets an extra pi/4 rotation to make it consistent with
% other carriers, as we need it for snr_update and scatter diagram
phase_difference(Nc+1) *= exp(j*pi/4);
else
% map (Nc,1) QPSK symbols back into an (1,Nc*Nb) array of bits
rx_bits(1:Nb:Nc*Nb) = real(rx_symbols) > 0;
rx_bits(2:Nb:Nc*Nb) = imag(rx_symbols) > 0;
endif
endfunction
% given phase differences update estimates of signal and noise levels
function [sig_est noise_est] = snr_update(sig_est, noise_est, phase_difference)
global snr_coeff;
global Nc;
% mag of each symbol is distance from origin, this gives us a
% vector of mags, one for each carrier.
s = abs(phase_difference);
% signal mag estimate for each carrier is a smoothed version
% of instantaneous magntitude, this gives us a vector of smoothed
% mag estimates, one for each carrier.
sig_est = snr_coeff*sig_est + (1 - snr_coeff)*s;
% noise mag estimate is distance of current symbol from average
% location of that symbol. We reflect all symbols into the first
% quadrant for convenience.
refl_symbols = abs(real(phase_difference)) + j*abs(imag(phase_difference));
n = abs(exp(j*pi/4)*sig_est - refl_symbols);
% noise mag estimate for each carrier is a smoothed version of
% instantaneous noise mag, this gives us a vector of smoothed
% noise power estimates, one for each carrier.
noise_est = snr_coeff*noise_est + (1 - snr_coeff)*n;
endfunction
% calculate current SNR estimate (3000Hz noise BW)
function snr_dB = calc_snr(sig_est, noise_est)
global Rs;
% find total signal power by summing power in all carriers
S = sum(sig_est .^2);
SdB = 10*log10(S);
% Average noise mag across all carriers and square to get an average
% noise power. This is an estimate of the noise power in Rs = 50Hz of
% BW (note for raised root cosine filters Rs is the noise BW of the
% filter)
N50 = mean(noise_est).^2;
N50dB = 10*log10(N50);
% Now multiply by (3000 Hz)/(50 Hz) to find the total noise power in
% 3000 Hz
N3000dB = N50dB + 10*log10(3000/Rs);
snr_dB = SdB - N3000dB;
endfunction
% returns nbits from a repeating sequence of random data
function bits = get_test_bits(nbits)
global Ntest_bits; % length of test sequence
global current_test_bit;
global test_bits;
for i=1:nbits
bits(i) = test_bits(current_test_bit++);
if (current_test_bit > Ntest_bits)
current_test_bit = 1;
endif
end
endfunction
% Accepts nbits from rx and attempts to sync with test_bits sequence.
% if sync OK measures bit errors
function [sync bit_errors] = put_test_bits(rx_bits)
global Ntest_bits; % length of test sequence
global test_bits;
global rx_test_bits_mem;
% Append to our memory
[m n] = size(rx_bits);
rx_test_bits_mem(1:Ntest_bits-n) = rx_test_bits_mem(n+1:Ntest_bits);
rx_test_bits_mem(Ntest_bits-n+1:Ntest_bits) = rx_bits;
% see how many bit errors we get when checked against test sequence
bit_errors = sum(xor(test_bits,rx_test_bits_mem));
% if less than a thresh we are aligned and in sync with test sequence
ber = bit_errors/Ntest_bits;
sync = 0;
if (ber < 0.2)
sync = 1;
endif
endfunction
% Generate M samples of DBPSK pilot signal for Freq offset estimation
function [pilot_fdm bit symbol filter_mem phase] = generate_pilot_fdm(bit, symbol, filter_mem, phase, freq)
global M;
global Nfilter;
global gt_alpha5_root;
% +1 -1 +1 -1 DBPSK sync carrier, once filtered becomes two spectral
% lines at +/- Rs/2
if bit
symbol = -symbol;
else
symbol = symbol;
end
if bit
bit = 0;
else
bit = 1;
end
% filter DPSK symbol to create M baseband samples
filter_mem(Nfilter) = (sqrt(2)/2)*symbol;
for i=1:M
tx_baseband(i) = M*filter_mem(M:M:Nfilter) * gt_alpha5_root(M-i+1:M:Nfilter)';
end
filter_mem(1:Nfilter-M) = filter_mem(M+1:Nfilter);
filter_mem(Nfilter-M+1:Nfilter) = zeros(1,M);
% upconvert
for i=1:M
phase = phase * freq;
pilot_fdm(i) = sqrt(2)*2*tx_baseband(i)*phase;
end
endfunction
% Generate a 4M sample vector of DBPSK pilot signal. As the pilot signal
% is periodic in 4M samples we can then use this vector as a look up table
% for pilot signal generation in the demod.
function pilot_lut = generate_pilot_lut()
global Nc;
global Nfilter;
global M;
global freq;
% pilot states
pilot_rx_bit = 0;
pilot_symbol = sqrt(2);
pilot_freq = freq(Nc+1);
pilot_phase = 1;
pilot_filter_mem = zeros(1, Nfilter);
%prev_pilot = zeros(M,1);
pilot_lut = [];
F=8;
for f=1:F
[pilot pilot_rx_bit pilot_symbol pilot_filter_mem pilot_phase] = generate_pilot_fdm(pilot_rx_bit, pilot_symbol, pilot_filter_mem, pilot_phase, pilot_freq);
%prev_pilot = pilot;
pilot_lut = [pilot_lut pilot];
end
% discard first 4 symbols as filter memory is filling, just keep last
% four symbols
pilot_lut = pilot_lut(4*M+1:M*F);
endfunction
% grab next pilot samples for freq offset estimation at demod
function [pilot prev_pilot pilot_lut_index prev_pilot_lut_index] = get_pilot(pilot_lut_index, prev_pilot_lut_index, nin)
global M;
global pilot_lut;
for i=1:nin
pilot(i) = pilot_lut(pilot_lut_index);
pilot_lut_index++;
if pilot_lut_index > 4*M
pilot_lut_index = 1;
end
prev_pilot(i) = pilot_lut(prev_pilot_lut_index);
prev_pilot_lut_index++;
if prev_pilot_lut_index > 4*M
prev_pilot_lut_index = 1;
end
end
endfunction
% Change the sample rate by a small amount, for example 1000ppm (ratio
% = 1.001). Always returns nout samples in buf_out, but uses a
% variable number of input samples nin to accomodate the change in
% sample rate. nin is nominally set to nout, but may use nout +/- 2
% samples to accomodate the different sample rates. buf_in should be
% of length nout+6 samples to accomodate this, and buf_in should be
% updated externally based on the nin returned each time. "ratio" is
% Fs_in/Fs_out, for example 48048/48000 = 1.001 (+1000ppm) or
% 47952/48000 = 0.999 (-1000ppm). Uses linear interpolation to
% perform the resampling. This requires a highly over-sampled signal,
% for example 48000Hz sample rate for the modem signal centred on
% 1kHz, otherwise linear interpolation will have a low pass filter effect
% (for example an 8000Hz sample rate for modem signal centred on 1kHz
% would cause problems).
function [buf_out t nin] = resample(buf_in, t, ratio, nout)
for i=1:nout
c = floor(t);
a = t - c;
b = 1 - a;
buf_out(i) = buf_in(c)*b + buf_in(c+1)*a;
t += ratio;
end
t -= nout;
% adjust nin and t so that on next call we start with 3 < t < 4,
% this gives us +/- 2 samples room to move before we hit start or
% end of buf_in
delta = floor(t - 3);
nin = nout + delta;
t -= delta;
endfunction
% freq offset state machine. Moves between acquire and track states based
% on BPSK pilot sequence. Freq offset estimator occasionally makes mistakes
% when used continuously. So we use it until we have acquired the BPSK pilot,
% then switch to a more robust tracking algorithm. If we lose sync we switch
% back to acquire mode for fast-requisition.
function [track state] = freq_state(sync_bit, state)
% acquire state, look for 6 symbol 010101 sequence from sync bit
next_state = state;
if state == 0
if sync_bit == 0
next_state = 1;
end
end
if state == 1
if sync_bit == 1
next_state = 2;
else
next_state = 0;
end
end
if state == 2
if sync_bit == 0
next_state = 3;
else
next_state = 0;
end
end
if state == 3
if sync_bit == 1
next_state = 4;
else
next_state = 0;
end
end
if state == 4
if sync_bit == 0
next_state = 5;
else
next_state = 0;
end
end
if state == 5
if sync_bit == 1
next_state = 6;
else
next_state = 0;
end
end
% states 6 and above are track mode, make sure we keep getting 0101 sync bit sequence
if state == 6
if sync_bit == 0
next_state = 7;
else
next_state = 0;
end
end
if state == 7
if sync_bit == 1
next_state = 6;
else
next_state = 0;
end
end
state = next_state;
if state >= 6
track = 1;
else
track = 0;
end
endfunction
% Save test bits to a text file in the form of a C array
function test_bits_file(filename)
global test_bits;
global Ntest_bits;
f=fopen(filename,"wt");
fprintf(f,"/* Generated by test_bits_file() Octave function */\n\n");
fprintf(f,"const int test_bits[]={\n");
for m=1:Ntest_bits-1
fprintf(f," %d,\n",test_bits(m));
endfor
fprintf(f," %d\n};\n",test_bits(Ntest_bits));
fclose(f);
endfunction
% Saves RN filter coeffs to a text file in the form of a C array
function rn_file(filename)
global gt_alpha5_root;
global Nfilter;
f=fopen(filename,"wt");
fprintf(f,"/* Generated by rn_file() Octave function */\n\n");
fprintf(f,"const float gt_alpha5_root[]={\n");
for m=1:Nfilter-1
fprintf(f," %g,\n",gt_alpha5_root(m));
endfor
fprintf(f," %g\n};\n",gt_alpha5_root(Nfilter));
fclose(f);
endfunction
function pilot_coeff_file(filename)
global pilot_coeff;
global Npilotcoeff;
f=fopen(filename,"wt");
fprintf(f,"/* Generated by pilot_coeff_file() Octave function */\n\n");
fprintf(f,"const float pilot_coeff[]={\n");
for m=1:Npilotcoeff-1
fprintf(f," %g,\n",pilot_coeff(m));
endfor
fprintf(f," %g\n};\n",pilot_coeff(Npilotcoeff));
fclose(f);
endfunction
% Saves hanning window coeffs to a text file in the form of a C array
function hanning_file(filename)
global Npilotlpf;
h = hanning(Npilotlpf);
f=fopen(filename,"wt");
fprintf(f,"/* Generated by hanning_file() Octave function */\n\n");
fprintf(f,"const float hanning[]={\n");
for m=1:Npilotlpf-1
fprintf(f," %g,\n", h(m));
endfor
fprintf(f," %g\n};\n", h(Npilotlpf));
fclose(f);
endfunction
function png_file(fig, pngfilename)
figure(fig);
pngname = sprintf("%s.png",pngfilename);
print(pngname, '-dpng', "-S500,500")
pngname = sprintf("%s_large.png",pngfilename);
print(pngname, '-dpng', "-S800,600")
endfunction
% Initialise ----------------------------------------------------
global pilot_bit;
pilot_bit = 0; % current value of pilot bit
global tx_filter_memory;
tx_filter_memory = zeros(Nc+1, Nfilter);
global rx_filter_memory;
rx_filter_memory = zeros(Nc+1, Nfilter);
% phasors used for up and down converters
global freq;
freq = zeros(Nc+1,1);
for c=1:Nc/2
carrier_freq = (-Nc/2 - 1 + c)*Fsep + Fcentre;
freq(c) = exp(j*2*pi*carrier_freq/Fs);
end
for c=Nc/2+1:Nc
carrier_freq = (-Nc/2 + c)*Fsep + Fcentre;
freq(c) = exp(j*2*pi*carrier_freq/Fs);
end
freq(Nc+1) = exp(j*2*pi*Fcentre/Fs);
% Spread initial FDM carrier phase out as far as possible. This
% helped PAPR for a few dB. We don't need to adjust rx phase as DQPSK
% takes care of that.
global phase_tx;
%phase_tx = ones(Nc+1,1);
phase_tx = exp(j*2*pi*(0:Nc)/(Nc+1));
global phase_rx;
phase_rx = ones(Nc+1,1);
% Freq offset estimator constants
global Mpilotfft = 256;
global Npilotcoeff = 30; % number of pilot LPF coeffs
global pilot_coeff = fir1(Npilotcoeff-1, 200/(Fs/2))'; % 200Hz LPF
global Npilotbaseband = Npilotcoeff + M + M/P; % number of pilot baseband samples reqd for pilot LPF
global Npilotlpf = 4*M; % number of samples we DFT pilot over, pilot est window
% pilot LUT, used for copy of pilot at rx
global pilot_lut;
pilot_lut = generate_pilot_lut();
pilot_lut_index = 1;
prev_pilot_lut_index = 3*M+1;
% Freq offset estimator states
global pilot_baseband1;
global pilot_baseband2;
pilot_baseband1 = zeros(1, Npilotbaseband); % pilot baseband samples
pilot_baseband2 = zeros(1, Npilotbaseband); % pilot baseband samples
global pilot_lpf1
global pilot_lpf2
pilot_lpf1 = zeros(1, Npilotlpf); % LPF pilot samples
pilot_lpf2 = zeros(1, Npilotlpf); % LPF pilot samples
% Timing estimator states
global rx_filter_mem_timing;
rx_filter_mem_timing = zeros(Nc+1, Nt*P);
global rx_baseband_mem_timing;
rx_baseband_mem_timing = zeros(Nc+1, Nfiltertiming);
% Test bit stream constants
global Ntest_bits = Nc*Nb*4; % length of test sequence
global test_bits = rand(1,Ntest_bits) > 0.5;
% Test bit stream state variables
global current_test_bit = 1;
current_test_bit = 1;
global rx_test_bits_mem;
rx_test_bits_mem = zeros(1,Ntest_bits);

217
libs/libcodec2/octave/fdmdv_demod.m Normal file
View File

@ -0,0 +1,217 @@
% fdmdv_demod.m
%
% Demodulator function for FDMDV modem (Octave version). Requires
% 8kHz sample rate raw files as input
%
% Copyright David Rowe 2012
% This program is distributed under the terms of the GNU General Public License
% Version 2
%
function fdmdv_demod(rawfilename, nbits, pngname)
fdmdv; % include modem code
modulation = 'dqpsk';
fin = fopen(rawfilename, "rb");
gain = 1000;
frames = nbits/(Nc*Nb);
prev_rx_symbols = ones(Nc+1,1);
foff_phase = 1;
% BER stats
total_bit_errors = 0;
total_bits = 0;
bit_errors_log = [];
sync_log = [];
test_frame_sync_log = [];
test_frame_sync_state = 0;
% SNR states
sig_est = zeros(Nc+1,1);
noise_est = zeros(Nc+1,1);
% logs of various states for plotting
rx_symbols_log = [];
rx_timing_log = [];
foff_log = [];
rx_fdm_log = [];
snr_est_log = [];
% misc states
nin = M; % timing correction for sample rate differences
foff = 0;
track_log = [];
track = 0;
fest_state = 0;
% Main loop ----------------------------------------------------
for f=1:frames
% obtain nin samples of the test input signal
for i=1:nin
rx_fdm(i) = fread(fin, 1, "short")/gain;
end
rx_fdm_log = [rx_fdm_log rx_fdm(1:nin)];
% frequency offset estimation and correction
[pilot prev_pilot pilot_lut_index prev_pilot_lut_index] = get_pilot(pilot_lut_index, prev_pilot_lut_index, nin);
[foff_coarse S1 S2] = rx_est_freq_offset(rx_fdm, pilot, prev_pilot, nin);
if track == 0
foff = foff_coarse;
end
foff_log = [ foff_log foff ];
foff_rect = exp(j*2*pi*foff/Fs);
for i=1:nin
foff_phase *= foff_rect';
rx_fdm(i) = rx_fdm(i)*foff_phase;
end
% baseband processing
rx_baseband = fdm_downconvert(rx_fdm, nin);
rx_filt = rx_filter(rx_baseband, nin);
[rx_symbols rx_timing] = rx_est_timing(rx_filt, rx_baseband, nin);
rx_timing_log = [rx_timing_log rx_timing];
nin = M;
if rx_timing > 2*M/P
nin += M/P;
end
if rx_timing < 0;
nin -= M/P;
end
if strcmp(modulation,'dqpsk')
rx_symbols_log = [rx_symbols_log rx_symbols.*conj(prev_rx_symbols)*exp(j*pi/4)];
else
rx_symbols_log = [rx_symbols_log rx_symbols];
endif
[rx_bits sync f_err pd] = qpsk_to_bits(prev_rx_symbols, rx_symbols, modulation);
[sig_est noise_est] = snr_update(sig_est, noise_est, pd);
snr_est = calc_snr(sig_est, noise_est);
snr_est_log = [snr_est_log snr_est];
foff -= 0.5*f_err;
prev_rx_symbols = rx_symbols;
sync_log = [sync_log sync];
% freq est state machine
[track fest_state] = freq_state(sync, fest_state);
track_log = [track_log track];
% count bit errors if we find a test frame
[test_frame_sync bit_errors] = put_test_bits(rx_bits);
if (test_frame_sync == 1)
total_bit_errors = total_bit_errors + bit_errors;
total_bits = total_bits + Ntest_bits;
bit_errors_log = [bit_errors_log bit_errors/Ntest_bits];
else
bit_errors_log = [bit_errors_log 0];
end
% test frame sync state machine, just for more informative plots
next_test_frame_sync_state = test_frame_sync_state;
if (test_frame_sync_state == 0)
if (test_frame_sync == 1)
next_test_frame_sync_state = 1;
test_frame_count = 0;
end
end
if (test_frame_sync_state == 1)
% we only expect another test_frame_sync pulse every 4 symbols
test_frame_count++;
if (test_frame_count == 4)
test_frame_count = 0;
if ((test_frame_sync == 0))
next_test_frame_sync_state = 0;
end
end
end
test_frame_sync_state = next_test_frame_sync_state;
test_frame_sync_log = [test_frame_sync_log test_frame_sync_state];
end
% ---------------------------------------------------------------------
% Print Stats
% ---------------------------------------------------------------------
ber = total_bit_errors / total_bits;
printf("%d bits %d errors BER: %1.4f\n",total_bits, total_bit_errors, ber);
% ---------------------------------------------------------------------
% Plots
% ---------------------------------------------------------------------
xt = (1:frames)/Rs;
secs = frames/Rs;
figure(1)
clf;
[n m] = size(rx_symbols_log);
plot(real(rx_symbols_log(1:Nc+1,15:m)),imag(rx_symbols_log(1:Nc+1,15:m)),'+')
axis([-2 2 -2 2]);
title('Scatter Diagram');
figure(2)
clf;
subplot(211)
plot(xt, rx_timing_log)
title('timing offset (samples)');
subplot(212)
plot(xt, foff_log, '-;freq offset;')
hold on;
plot(xt, track_log*75, 'r;course-fine;');
hold off;
title('Freq offset (Hz)');
grid
figure(3)
clf;
subplot(211)
[a b] = size(rx_fdm_log);
xt1 = (1:b)/Fs;
plot(xt1, rx_fdm_log);
title('Rx FDM Signal');
subplot(212)
spec(rx_fdm_log,8000);
title('FDM Rx Spectrogram');
figure(4)
clf;
subplot(311)
stem(xt, sync_log)
axis([0 secs 0 1.5]);
title('BPSK Sync')
subplot(312)
stem(xt, bit_errors_log);
title('Bit Errors for test frames')
subplot(313)
plot(xt, test_frame_sync_log);
axis([0 secs 0 1.5]);
title('Test Frame Sync')
figure(5)
clf;
plot(xt, snr_est_log);
title('SNR Estimates')
endfunction

128
libs/libcodec2/octave/fdmdv_demod_c.m Normal file
View File

@ -0,0 +1,128 @@
% fdmdv_demod_c.m
%
% Plots Octave dump file information from C FDMDV demodulator program,
% to give a similar set of plots to fdmdv_demod.m. Useful for off
% line analysis of demod performance.
%
% Copyright David Rowe 2012
% This program is distributed under the terms of the GNU General Public License
% Version 2
%
function fdmdv_demod_c(dumpfilename, bits)
fdmdv; % include modem code
frames = bits/(Nc*Nb);
load(dumpfilename);
% BER stats
total_bit_errors = 0;
total_bits = 0;
bit_errors_log = [];
sync_log = [];
test_frame_sync_log = [];
test_frame_sync_state = 0;
% Run thru received bits to look for test pattern
bits_per_frame = Nc*Nb;
for f=1:frames
rx_bits = rx_bits_log_c((f-1)*bits_per_frame+1:f*bits_per_frame);
% count bit errors if we find a test frame
[test_frame_sync bit_errors] = put_test_bits(rx_bits);
if (test_frame_sync == 1)
total_bit_errors = total_bit_errors + bit_errors;
total_bits = total_bits + Ntest_bits;
bit_errors_log = [bit_errors_log bit_errors/Ntest_bits];
else
bit_errors_log = [bit_errors_log 0];
end
% test frame sync state machine, just for more informative plots
next_test_frame_sync_state = test_frame_sync_state;
if (test_frame_sync_state == 0)
if (test_frame_sync == 1)
next_test_frame_sync_state = 1;
test_frame_count = 0;
end
end
if (test_frame_sync_state == 1)
% we only expect another test_frame_sync pulse every 4 symbols
test_frame_count++;
if (test_frame_count == 4)
test_frame_count = 0;
if ((test_frame_sync == 0))
next_test_frame_sync_state = 0;
end
end
end
test_frame_sync_state = next_test_frame_sync_state;
test_frame_sync_log = [test_frame_sync_log test_frame_sync_state];
end
% ---------------------------------------------------------------------
% Plots
% ---------------------------------------------------------------------
xt = (1:frames)/Rs;
secs = frames/Rs;
figure(1)
clf;
plot(real(rx_symbols_log_c(1:Nc+1,15:frames)),imag(rx_symbols_log_c(1:Nc+1,15:frames)),'+')
axis([-2 2 -2 2]);
title('Scatter Diagram');
figure(2)
clf;
subplot(211)
plot(xt, rx_timing_log_c(1:frames))
title('timing offset (samples)');
subplot(212)
plot(xt, foff_log_c(1:frames), '-;freq offset;')
hold on;
plot(xt, coarse_fine_log_c(1:frames)*75, 'r;course-fine;');
hold off;
title('Freq offset (Hz)');
grid
figure(3)
clf;
subplot(211)
b = M*frames;
xt1 = (1:b)/Fs;
plot(xt1, rx_fdm_log_c(1:b));
title('Rx FDM Signal');
subplot(212)
spec(rx_fdm_log_c(1:b),8000);
title('FDM Rx Spectrogram');
figure(4)
clf;
subplot(311)
stem(xt, sync_bit_log_c(1:frames))
axis([0 secs 0 1.5]);
title('BPSK Sync')
subplot(312)
stem(xt, bit_errors_log);
title('Bit Errors for test frames')
subplot(313)
plot(xt, test_frame_sync_log);
axis([0 secs 0 1.5]);
title('Test Frame Sync')
figure(5)
clf;
plot(xt, snr_est_log_c(1:frames));
title('SNR Estimates')
endfunction

32
libs/libcodec2/octave/fdmdv_mod.m Normal file
View File

@ -0,0 +1,32 @@
% fdmdv_mod.m
%
% Modulator function for FDMDV modem, uses test frames as input and
% outputs a raw file of 16 bit shorts at a sample rate of 8 kHz.
%
% Copyright David Rowe 2012
% This program is distributed under the terms of the GNU General Public License
% Version 2
%
function tx_fdm = fdmdv_mod(rawfilename, nbits)
fdmdv; % include modem code
frames = floor(nbits/(Nc*Nb))
tx_fdm = [];
gain = 1000; % Scale up to 16 bit shorts
prev_tx_symbols = ones(Nc+1,1);
for i=1:frames
tx_bits = get_test_bits(Nc*Nb);
tx_symbols = bits_to_qpsk(prev_tx_symbols, tx_bits,'dqpsk');
prev_tx_symbols = tx_symbols;
tx_baseband = tx_filter(tx_symbols);
tx_fdm = [tx_fdm real(fdm_upconvert(tx_baseband))];
end
tx_fdm *= gain;
fout = fopen(rawfilename,"wb");
fwrite(fout, tx_fdm, "short");
fclose(fout);
endfunction

318
libs/libcodec2/octave/fdmdv_ut.m Normal file
View File

@ -0,0 +1,318 @@
% fdmdv_ut.m
%
% Unit Test program for FDMDV modem. Useful for general development as it has
% both tx and rx sides, and basic AWGN channel simulation.
%
% Copyright David Rowe 2012
% This program is distributed under the terms of the GNU General Public License
% Version 2
%
fdmdv; % load modem code
% Simulation Parameters --------------------------------------
frames = 25;
EbNo_dB = 7.3;
Foff_hz = 0;
modulation = 'dqpsk';
hpa_clip = 150;
% ------------------------------------------------------------
tx_filt = zeros(Nc,M);
rx_symbols_log = [];
rx_phase_log = 0;
rx_timing_log = 0;
tx_pwr = 0;
noise_pwr = 0;
rx_fdm_log = [];
rx_baseband_log = [];
rx_bits_offset = zeros(Nc*Nb*2);
prev_tx_symbols = ones(Nc+1,1);
prev_rx_symbols = ones(Nc+1,1);
ferr = 0;
foff = 0;
foff_log = [];
tx_baseband_log = [];
tx_fdm_log = [];
% BER stats
total_bit_errors = 0;
total_bits = 0;
bit_errors_log = [];
sync_log = [];
test_frame_sync_log = [];
test_frame_sync_state = 0;
% SNR estimation states
sig_est = zeros(Nc+1,1);
noise_est = zeros(Nc+1,1);
% fixed delay simuation
Ndelay = M+20;
rx_fdm_delay = zeros(Ndelay,1);
% ---------------------------------------------------------------------
% Eb/No calculations. We need to work out Eb/No for each FDM carrier.
% Total power is sum of power in all FDM carriers
% ---------------------------------------------------------------------
C = 1; % power of each FDM carrier (energy/sample). Total Carrier power should = Nc*C = Nc
N = 1; % total noise power (energy/sample) of noise source across entire bandwidth
% Eb = Carrier power * symbol time / (bits/symbol)
% = C *(1/Rs) / 2
Eb_dB = 10*log10(C) - 10*log10(Rs) - 10*log10(2);
No_dBHz = Eb_dB - EbNo_dB;
% Noise power = Noise spectral density * bandwidth
% Noise power = Noise spectral density * Fs/2 for real signals
N_dB = No_dBHz + 10*log10(Fs/2);
Ngain_dB = N_dB - 10*log10(N);
Ngain = 10^(Ngain_dB/20);
% C/No = Carrier Power/noise spectral density
% = power per carrier*number of carriers / noise spectral density
CNo_dB = 10*log10(C) + 10*log10(Nc) - No_dBHz;
% SNR in equivalent 3000 Hz SSB channel
B = 3000;
SNR = CNo_dB - 10*log10(B);
% freq offset simulation states
phase_offset = 1;
freq_offset = exp(j*2*pi*Foff_hz/Fs);
foff_phase = 1;
t = 0;
foff = 0;
fest_state = 0;
track = 0;
track_log = [];
% ---------------------------------------------------------------------
% Main loop
% ---------------------------------------------------------------------
for f=1:frames
% -------------------
% Modulator
% -------------------
tx_bits = get_test_bits(Nc*Nb);
tx_symbols = bits_to_qpsk(prev_tx_symbols, tx_bits, modulation);
prev_tx_symbols = tx_symbols;
tx_baseband = tx_filter(tx_symbols);
tx_baseband_log = [tx_baseband_log tx_baseband];
tx_fdm = fdm_upconvert(tx_baseband);
tx_pwr = 0.9*tx_pwr + 0.1*real(tx_fdm)*real(tx_fdm)'/(M);
% -------------------
% Channel simulation
% -------------------
% frequency offset
%Foff_hz += 1/Rs;
Foff = Foff_hz;
for i=1:M
% Time varying freq offset
%Foff = Foff_hz + 100*sin(t*2*pi/(300*Fs));
%t++;
freq_offset = exp(j*2*pi*Foff/Fs);
phase_offset *= freq_offset;
tx_fdm(i) = phase_offset*tx_fdm(i);
end
tx_fdm = real(tx_fdm);
% HPA non-linearity
tx_fdm(find(abs(tx_fdm) > hpa_clip)) = hpa_clip;
tx_fdm_log = [tx_fdm_log tx_fdm];
rx_fdm = tx_fdm;
% AWGN noise
noise = Ngain*randn(1,M);
noise_pwr = 0.9*noise_pwr + 0.1*noise*noise'/M;
rx_fdm += noise;
rx_fdm_log = [rx_fdm_log rx_fdm];
% Delay
rx_fdm_delay(1:Ndelay-M) = rx_fdm_delay(M+1:Ndelay);
rx_fdm_delay(Ndelay-M+1:Ndelay) = rx_fdm;
%rx_fdm_delay = rx_fdm;
% -------------------
% Demodulator
% -------------------
% frequency offset estimation and correction, need to call rx_est_freq_offset even in track
% mode to keep states updated
[pilot prev_pilot pilot_lut_index prev_pilot_lut_index] = get_pilot(pilot_lut_index, prev_pilot_lut_index, M);
[foff_course S1 S2] = rx_est_freq_offset(rx_fdm_delay, pilot, prev_pilot, M);
if track == 0
foff = foff_course;
end
foff_log = [ foff_log foff ];
foff_rect = exp(j*2*pi*foff/Fs);
for i=1:M
foff_phase *= foff_rect';
rx_fdm_delay(i) = rx_fdm_delay(i)*foff_phase;
end
% baseband processing
rx_baseband = fdm_downconvert(rx_fdm_delay(1:M), M);
rx_baseband_log = [rx_baseband_log rx_baseband];
rx_filt = rx_filter(rx_baseband, M);
[rx_symbols rx_timing] = rx_est_timing(rx_filt, rx_baseband, M);
rx_timing_log = [rx_timing_log rx_timing];
%rx_phase = rx_est_phase(rx_symbols);
%rx_phase_log = [rx_phase_log rx_phase];
%rx_symbols = rx_symbols*exp(j*rx_phase);
[rx_bits sync foff_fine pd] = qpsk_to_bits(prev_rx_symbols, rx_symbols, modulation);
if strcmp(modulation,'dqpsk')
%rx_symbols_log = [rx_symbols_log rx_symbols.*conj(prev_rx_symbols)*exp(j*pi/4)];
rx_symbols_log = [rx_symbols_log pd];
else
rx_symbols_log = [rx_symbols_log rx_symbols];
endif
foff -= 0.5*ferr;
prev_rx_symbols = rx_symbols;
sync_log = [sync_log sync];
% freq est state machine
[track fest_state] = freq_state(sync, fest_state);
track_log = [track_log track];
% Update SNR est
[sig_est noise_est] = snr_update(sig_est, noise_est, pd);
% count bit errors if we find a test frame
% Allow 15 frames for filter memories to fill and time est to settle
[test_frame_sync bit_errors] = put_test_bits(rx_bits);
if test_frame_sync == 1
total_bit_errors = total_bit_errors + bit_errors;
total_bits = total_bits + Ntest_bits;
bit_errors_log = [bit_errors_log bit_errors];
else
bit_errors_log = [bit_errors_log 0];
end
% test frame sync state machine, just for more informative plots
next_test_frame_sync_state = test_frame_sync_state;
if (test_frame_sync_state == 0)
if (test_frame_sync == 1)
next_test_frame_sync_state = 1;
test_frame_count = 0;
end
end
if (test_frame_sync_state == 1)
% we only expect another test_frame_sync pulse every 4 symbols
test_frame_count++;
if (test_frame_count == 4)
test_frame_count = 0;
if ((test_frame_sync == 0))
next_test_frame_sync_state = 0;
end
end
end
test_frame_sync_state = next_test_frame_sync_state;
test_frame_sync_log = [test_frame_sync_log test_frame_sync_state];
end
% ---------------------------------------------------------------------
% Print Stats
% ---------------------------------------------------------------------
ber = total_bit_errors / total_bits;
% Peak to Average Power Ratio calcs from http://www.dsplog.com
papr = max(tx_fdm_log.*conj(tx_fdm_log)) / mean(tx_fdm_log.*conj(tx_fdm_log));
papr_dB = 10*log10(papr);
% Note Eb/No set point is for Nc data carriers only, exclduing pilot.
% This is convenient for testing BER versus Eb/No. Measured Eb/No
% includes power of pilot. Similar for SNR, first number is SNR excluding
% pilot pwr for Eb/No set point, 2nd value is measured SNR which will be a little
% higher as pilot power is included.
printf("Eb/No (meas): %2.2f (%2.2f) dB\n", EbNo_dB, 10*log10(0.25*tx_pwr*Fs/(Rs*Nc*noise_pwr)));
printf("bits........: %d\n", total_bits);
printf("errors......: %d\n", total_bit_errors);
printf("BER.........: %1.4f\n", ber);
printf("PAPR........: %1.2f dB\n", papr_dB);
printf("SNR...(meas): %2.2f (%2.2f) dB\n", SNR, calc_snr(sig_est, noise_est));
% ---------------------------------------------------------------------
% Plots
% ---------------------------------------------------------------------
figure(1)
clf;
[n m] = size(rx_symbols_log);
plot(real(rx_symbols_log(1:Nc+1,15:m)),imag(rx_symbols_log(1:Nc+1,15:m)),'+')
axis([-3 3 -3 3]);
title('Scatter Diagram');
figure(2)
clf;
subplot(211)
plot(rx_timing_log)
title('timing offset (samples)');
subplot(212)
plot(foff_log, '-;freq offset;')
hold on;
plot(track_log*75, 'r;course-fine;');
hold off;
title('Freq offset (Hz)');
figure(3)
clf;
subplot(211)
plot(real(tx_fdm_log));
title('FDM Tx Signal');
subplot(212)
Nfft=Fs;
S=fft(rx_fdm_log,Nfft);
SdB=20*log10(abs(S));
plot(SdB(1:Fs/4))
title('FDM Rx Spectrum');
figure(4)
clf;
subplot(311)
stem(sync_log)
axis([0 frames 0 1.5]);
title('BPSK Sync')
subplot(312)
stem(bit_errors_log);
title('Bit Errors for test frames')
subplot(313)
plot(test_frame_sync_log);
axis([0 frames 0 1.5]);
title('Test Frame Sync')

40
libs/libcodec2/octave/gen_rn_coeffs.m Normal file
View File

@ -0,0 +1,40 @@
% gen_rn_coeffs.m
% David Rowe 13 april 2012
%
% Generate root raised cosine (Root Nyquist) filter coefficients
% thanks http://www.dsplog.com/db-install/wp-content/uploads/2008/05/raised_cosine_filter.m
function coeffs = gen_rn_coeffs(alpha, T, Rs, Nsym, M)
Ts = 1/Rs;
n = -Nsym*Ts/2:T:Nsym*Ts/2;
Nfilter = Nsym*M;
Nfiltertiming = M+Nfilter+M;
sincNum = sin(pi*n/Ts); % numerator of the sinc function
sincDen = (pi*n/Ts); % denominator of the sinc function
sincDenZero = find(abs(sincDen) < 10^-10);
sincOp = sincNum./sincDen;
sincOp(sincDenZero) = 1; % sin(pix/(pix) =1 for x =0
cosNum = cos(alpha*pi*n/Ts);
cosDen = (1-(2*alpha*n/Ts).^2);
cosDenZero = find(abs(cosDen)<10^-10);
cosOp = cosNum./cosDen;
cosOp(cosDenZero) = pi/4;
gt_alpha5 = sincOp.*cosOp;
Nfft = 4096;
GF_alpha5 = fft(gt_alpha5,Nfft)/M;
% sqrt causes stop band to be amplified, this hack pushes it down again
for i=1:Nfft
if (abs(GF_alpha5(i)) < 0.02)
GF_alpha5(i) *= 0.001;
endif
end
GF_alpha5_root = sqrt(abs(GF_alpha5)) .* exp(j*angle(GF_alpha5));
ifft_GF_alpha5_root = ifft(GF_alpha5_root);
coeffs = real((ifft_GF_alpha5_root(1:Nfilter)));
endfunction

12
libs/libcodec2/octave/hp_filt.m Normal file
View File

@ -0,0 +1,12 @@
% hp_filt.m
% David Rowe 20 Feb 2012
function hp_filt(in_file, out_file)
fin = fopen(in_file,"rb");
s = fread(fin,Inf,"short");
b = fir1(256, 300/4000, "high");
freqz(b);
s_hpf = filter(b,1,s);
fout = fopen(out_file,"wb");
fwrite(fout, s_hpf, "short");
endfunction

46
libs/libcodec2/octave/lpcpf.m Normal file
View File

@ -0,0 +1,46 @@
% lpcpf.m
% David Rowe Aug 27 2012
% Experiments with LPC post filtering
function lpcpf(ak_filename, f)
aks = load(ak_filename);
ak = aks(f,:);
[tmp p] = size(ak);
p -= 1;
A = freqz(1,ak, 4000);
AdB = 20*log10(abs(A));
gamma = 0.5;
gammas = gamma .^ (0:p);
W = freqz(ak .* gammas,1, 4000);
WdB = 20*log10(abs(W));
beta = 0.2;
R = abs(freqz(ak .* gammas, ak, 4000));
%P = (R/max(R)) .^ beta;
P = R .^ beta;
AP = abs(A) .* P;
eA = sum(abs(A) .^ 2);
eAP = sum(AP .^ 2);
gain = sqrt(eA/eAP)
AP *= gain;
PdB = 20*log10(P);
APdB = 20*log10(AP);
10*log10(sum(AP .^ 2))/10*log10(sum(abs(A) .^ 2))
figure(1);
clf;
plot(AdB);
hold on;
plot(WdB,'g');
plot(PdB,'r');
plot(APdB,'b.');
hold off;
endfunction

40
libs/libcodec2/octave/lspwarp.m Normal file
View File

@ -0,0 +1,40 @@
% lspwarp.m
% David Rowe Sep 2012
%
% Experimenting with non-linear LSP frequency axis for LSP quantisation
% Plots a scaled mel axis.
1;
function mel = freq2mel(f)
mel = 70*log10(1 + f/700);
endfunction
function freq = mel2freq(m)
freq = 700*(10 ^ (m/70) - 1);
endfunction
x = []; y = [];
for freq = 100:25:4000
mel = freq2mel(freq);
x = [x freq];
y = [y mel];
end
plot(x,y)
grid
mel_start = floor(freq2mel(100));
mel_end = floor(freq2mel(4000));
x = []; y = [];
for mel=mel_start:mel_end
freq = mel2freq(mel);
x = [x freq];
y = [y mel];
end
hold on;
plot(x,y, '+')
hold off;

47
libs/libcodec2/octave/phasesecord.m Normal file
View File

@ -0,0 +1,47 @@
% phasesecord.m
% David Rowe Aug 2012
% Used to experiment with aproximations of phase of 2nd order systems
function phasesecord(w,beta)
a = [1 -2*cos(w)*beta beta*beta];
b = 1;
[h w1] = freqz(b,a);
figure(1)
subplot(211)
plot(abs(h))
subplot(212)
plot(angle(h))
% for beta close to 1, we approximate 3 dB points as 1-beta above
% and below the resonance freq. Note this fails if w=0 as there is a
% double pole. Lets sample the freq response at the 3dB points and
% w:
ws = [w-(1-beta) w w+(1-beta)];
[h w1] = freqz(b,a,ws);
% gain as a fraction of max, should be 3dB. Within 1.3 dB or for w > pi/8,
% gets innacurate near w=0 due to 2nd pole
printf("mag measured...:"); printf("% 4.3f ", abs(h)/max(abs(h)));
% measured angle, 45 deg from angle at w
printf("\nangle measured.: "); printf("% 5.3f ", angle(h));
% Our estimate of angle, (pi+w) is phase at resonance, at lower 3dB
% phase is pi/4 ahead, at upper 3B pi/4 behind. -pi/2 is contribution of
% other pole at at -w to phase
ph_lower = (pi+w) + pi/4 - pi/2;
ph_res =(pi+w) - pi/2;
ph_upper = (pi+w) - pi/4 - pi/2;
ph_ests = [ph_lower ph_res ph_upper];
ph_ests = ph_ests - 2*pi*(floor(ph_ests/(2*pi)) + 0.5);
printf("\nangle estimated:"); printf("% 5.3f ", ph_ests);
printf("\n");
endfunction

150
libs/libcodec2/octave/pllpcpf.m Normal file
View File

@ -0,0 +1,150 @@
% Copyright David Rowe 2012
% This program is distributed under the terms of the GNU General Public License
% Version 2
%
% Plot amplitude modelling information from dump files to test and develop
% LPC post filter.
function pllpcpf(samname, f)
% switch some stuff off to unclutter display
plot_Am = 0;
plot_Amq = 0;
plot_err = 0;
plot_lsp = 0;
plot_snr = 0;
plot_vsnr = 0;
plot_sw = 0;
plot_pw = 1;
plot_pwb = 1;
plot_rw = 1;
sn_name = strcat(samname,"_sn.txt");
Sn = load(sn_name);
sw_name = strcat(samname,"_sw.txt");
Sw = load(sw_name);
sw__name = strcat(samname,"_sw_.txt");
if (file_in_path(".",sw__name))
Sw_ = load(sw__name);
endif
model_name = strcat(samname,"_model.txt");
model = load(model_name);
modelq_name = strcat(samname,"_qmodel.txt");
if (file_in_path(".",modelq_name))
modelq = load(modelq_name);
endif
% Pw (LPC synth filter spectrum) before post filter
pwb_name = strcat(samname,"_pwb.txt");
if (file_in_path(".",pwb_name))
Pwb = load(pwb_name);
endif
% Rw (Post filter spectrum)
rw_name = strcat(samname,"_rw.txt");
if (file_in_path(".",rw_name))
Rw = load(rw_name);
endif
% Pw (LPC synth filter spectrum) after post filter
pw_name = strcat(samname,"_pw.txt");
if (file_in_path(".",pw_name))
Pw = load(pw_name);
endif
Ew_on = 1;
k = ' ';
do
figure(1);
clf;
s = [ Sn(2*f-1,:) Sn(2*f,:) ];
size(s);
plot(s);
axis([1 length(s) -20000 20000]);
figure(2);
clf;
Wo = model(f,1);
L = model(f,2);
Am = model(f,3:(L+2));
if plot_Am
plot((1:L)*Wo*4000/pi, 20*log10(Am),";Am;r");
end
axis([1 4000 -10 80]);
hold on;
if plot_sw
plot((0:255)*4000/256, Sw(f,:),";Sw;");
end
if (file_in_path(".",modelq_name))
Amq = modelq(f,3:(L+2));
if plot_Amq
plot((1:L)*Wo*4000/pi, 20*log10(Amq),";Amq;g" );
end
if (file_in_path(".",pwb_name) && plot_pwb)
plot((0:255)*4000/256, 10*log10(Pwb(f,:)),";Pwb;r");
endif
if (file_in_path(".",rw_name) && plot_rw)
plot((0:255)*4000/256, 10*log10(Rw(f,:)),";Rw;b");
endif
if (file_in_path(".",pw_name) && plot_pw)
plot((0:255)*4000/256, 10*log10(Pw(f,:)),";Pw;g.");
endif
signal = Am * Am';
noise = (Am-Amq) * (Am-Amq)';
snr1 = 10*log10(signal/noise);
Am_err_label = sprintf(";Am error SNR %4.2f dB;m",snr1);
if plot_err
plot((1:L)*Wo*4000/pi, 20*log10(Amq) - 20*log10(Am), Am_err_label);
end
endif
hold off;
% interactive menu
printf("\rframe: %d menu: n-next b-back p-png q-quit", f);
fflush(stdout);
k = kbhit();
if (k == 'n')
f = f + 1;
endif
if (k == 'b')
f = f - 1;
endif
% optional print to PNG
if (k == 'p')
figure(1);
pngname = sprintf("%s_%d_sn.png",samname,f);
print(pngname, '-dpng', "-S500,500")
pngname = sprintf("%s_%d_sn_large.png",samname,f);
print(pngname, '-dpng', "-S800,600")
figure(2);
pngname = sprintf("%s_%d_sw.png",samname,f);
print(pngname, '-dpng', "-S500,500")
pngname = sprintf("%s_%d_sw_large.png",samname,f);
print(pngname, '-dpng', "-S1200,800")
endif
until (k == 'q')
printf("\n");
endfunction

46
libs/libcodec2/octave/pllsp.m Normal file
View File

@ -0,0 +1,46 @@
% Copyright David Rowe 2010
% This program is distributed under the terms of the GNU General Public License
% Version 2
%
% Plots a bunch of information related to LSP quantisation:
% - speech file
% - LSPs before and after quantisation
% - SNR for each frame
%
% Note: there is a 160 sample (two frame delay) from the when a sample
% enters the input buffer until it is at the centre of the analysis window
function pllsp(rawfile,
dumpfile_prefix_lpc_only,
dumpfile_prefix_lsp,
start_f, end_f)
fs=fopen(rawfile,"rb");
s=fread(fs,Inf,"short");
lpc_snr_name = strcat(dumpfile_prefix_lpc_only,"_lpc_snr.txt");
lpc10_snr = load(lpc_snr_name);
lpc_snr_name = strcat(dumpfile_prefix_lsp,"_lpc_snr.txt");
lsp_snr = load(lpc_snr_name);
lsp_name = strcat(dumpfile_prefix_lsp,"_lsp.txt");
lsps = load(lsp_name);
[m,n]=size(lsps);
lsp = lsps(1:2:m,:);
lsp_ = lsps(2:2:m,:);
figure(1);
clf;
subplot(211);
sp = s((start_f-2)*80:(end_f-2)*80);
plot(sp);
subplot(212);
plot(lpc10_snr((start_f+1):end_f)-lsp_snr((start_f+1):end_f));
figure(2);
plot((4000/pi)*lsp((start_f+1):end_f,:));
hold on;
plot((4000/pi)*lsp_((start_f+1):end_f,:),'+-');
hold off;
endfunction

27
libs/libcodec2/octave/pllspdt.m Normal file
View File

@ -0,0 +1,27 @@
% pllspdt.m
% Copyright David Rowe 2010
% This program is distributed under the terms of the GNU General Public License
% Version 2
%
% Test script to plot differences in LSps between frames
function pllspdt(rawfile,dumpfile_prefix_lsp,lspn, start_f, end_f)
fs=fopen(rawfile,"rb");
s=fread(fs,Inf,"short");
lsp_name = strcat(dumpfile_prefix_lsp,"_lsp.txt");
lsps = load(lsp_name);
[m,n]=size(lsps);
lsp = lsps(1:2:m,:);
lsp_ = lsps(2:2:m,:);
lspdt = lsp(2:m/2,:) - lsp(1:m/2-1,:);
figure(1);
clf;
sp = s((start_f-2)*80:(end_f-2)*80);
plot(sp);
figure(2);
plot((4000/pi)*lspdt((start_f+1):end_f,lspn));
endfunction

65
libs/libcodec2/octave/plppe.m Normal file
View File

@ -0,0 +1,65 @@
% Copyright David Rowe 2009
% This program is distributed under the terms of the GNU General Public License
% Version 2
%
% Plot two sparse phase prediction error text files.
% Generate data from print_pred_error, print_pred_error_sparse_wo_correction1 etc
function plppe(ppe1_file, ppe2_file, f)
ppe1 = load(ppe1_file);
ppe2 = load(ppe2_file);
std1 = std(nonzeros(ppe1(:,40:80)));
std2 = std(nonzeros(ppe2(:,40:80)));
printf("std dev for %s is %4.3f\n", ppe1_file, std1);
printf("std dev for %s is %4.3f\n", ppe2_file, std2);
figure(1);
clf;
subplot(211)
hist(nonzeros(ppe1(:,40:80)),20);
subplot(212)
hist(nonzeros(ppe2(:,40:80)),20);
k = ' ';
do
figure(2);
clf;
subplot(211)
L = length(nonzeros(ppe1(f,:)));
x = (1:L)*4000/L;
std1 = std(nonzeros(ppe1(f,:)));
legend = sprintf(";std dev %4.3f;", std1);
plot(x, nonzeros(ppe1(f,:)),legend);
axis([0 4000 -pi pi]);
subplot(212)
std2 = std(nonzeros(ppe2(f,:)));
legend = sprintf(";std dev %4.3f;", std2);
plot(x, nonzeros(ppe2(f,:)),legend);
axis([0 4000 -pi pi]);
% interactive menu
printf("\rframe: %d menu: n-next b-back p-png q-quit ", f);
fflush(stdout);
k = kbhit();
if (k == 'n')
f = f + 1;
endif
if (k == 'b')
f = f - 1;
endif
% optional print to PNG
if (k == 'p')
pngname = sprintf("%s_%d",samname,f);
png(pngname);
endif
until (k == 'q')
printf("\n");
endfunction

35
libs/libcodec2/octave/plsub.m Normal file
View File

@ -0,0 +1,35 @@
% Copyright David Rowe 2010
% This program is distributed under the terms of the GNU General Public License
% Version 2
% plots the difference of two files
function plsub(samname1, samname2, start_sam, end_sam, pngname)
fs1=fopen(samname1,"rb");
s1=fread(fs1,Inf,"short");
fs2=fopen(samname2,"rb");
s2=fread(fs2,Inf,"short");
st = 1;
en = length(s1);
if (nargin >= 3)
st = start_sam;
endif
if (nargin >= 4)
en = end_sam;
endif
figure(1);
clf;
l1 = strcat("r;",samname1,";");
plot(s1(st:en) - s2(st:en), l1);
%axis([1 en-st min(s1(st:en)) max(s1(st:en))]);
if (nargin == 5)
pngname = sprintf("%s.png",pngname);
print(pngname, '-dpng', "-S500,500")
pngname = sprintf("%s_large.png",pngname);
print(pngname, '-dpng', "-S800,600")
endif
endfunction

89
libs/libcodec2/octave/plvoicing.m Normal file
View File

@ -0,0 +1,89 @@
% Copyright David Rowe 2009
% This program is distributed under the terms of the GNU General Public License
% Version 2
%
% Plot voicing information from sample and dump files.
%
% samfilename is the raw source file, e.g. "../raw/hts1a.raw"
% samname is the dumpfile prefix, e.g. "../src/hts1a"
%
% There is a 160 sample (two frame delay) from the when a sample
% enters the input buffer until it is at the centre of the analysis window
function plvoicing(samfilename, samname, start_f, end_f, pngname)
fs=fopen(samfilename,"rb");
s=fread(fs,Inf,"short");
snr_name = strcat(samname,"_snr.txt");
snr = load(snr_name);
model_name = strcat(samname,"_model.txt");
model = load(model_name);
Wo = model((start_f+1):end_f,1);
F0 = Wo*4000/pi;
dF0 = F0(1:length(Wo)-1) - F0(2:length(Wo));
% work out LP and HP energy
for f=(start_f+1):end_f
L = model(f,2);
Am = model(f,3:(L+2));
L2 = floor(L/2);
elow = Am(1:L2) * Am(1:L2)';
ehigh = Am(L2:L) * Am(L2:L)';
erat(f-(start_f+1)+1) = 10*log10(elow/ehigh);
endfor
figure(1);
clf;
sp = s((start_f-2)*80:(end_f-2)*80);
plot(sp);
hold on;
vhigh = snr((start_f+1):end_f) > 7;
vlow = snr((start_f+1):end_f) > 4;
% test correction based on erat
vlowadj = vlow;
for f=1:length(erat)-1
if (vlow(f) == 0)
if (erat(f) > 10)
vlowadj(f) = 1;
endif
endif
if (vlow(f) == 1)
if (erat(f) < -10)
vlowadj(f) = 0;
endif
if (abs(dF0(f)) > 15)
vlowadj(f) = 0;
endif
endif
endfor
x = 1:(end_f-start_f);
plot(x*80,snr((start_f+1):end_f)*1000,';SNRdB x 1000;g+');
plot(x*80,-8000 + vhigh*2000,';7dB thresh;g');
plot(x*80,-11000 + vlowadj*2000,';vlow with corr;g');
plot(x*80,erat*1000,';elow/ehigh in dB;r');
plot(x*80,-14000 + vlow*2000,';4dB thresh;r');
hold off;
grid
if (nargin == 5)
print(pngname, "-dpng", "-S500,500")
endif
figure(2)
Wo = model((start_f+1):end_f,1);
F0 = Wo*4000/pi;
dF0 = F0(1:length(Wo)-1) - F0(2:length(Wo));
%plot(dF0,'+--')
%hold on;
%plot([ 1 length(dF0) ], [10 10] ,'r')
%plot([ 1 length(dF0) ], [-10 -10] ,'r')
%axis([1 length(dF0) -50 50])
%hold off;
plot(F0,'+--')
endfunction

101
libs/libcodec2/octave/sd.m Normal file
View File

@ -0,0 +1,101 @@
% sd.m
% David Rowe Aug 2012
% Plots the spectal distorion between twofiles of LPCs. Used for LSP
% quantisation tuning.
function sd(raw_filename, dump_file_prefix, f)
ak1_filename = sprintf("%s_ak.txt", dump_file_prefix);
ak2_filename = sprintf("%s_ak_.txt", dump_file_prefix);
ak1 = load(ak1_filename);
ak2 = load(ak2_filename);
[ak1_r, ak1_c] = size(ak1);
[ak2_r, ak2_c] = size(ak1);
frames = max([ak1_r ak2_r]);
sd = zeros(1,frames);
Ndft = 512;
A1 = zeros(frames, Ndft);
A2 = zeros(frames, Ndft);
% initial helicopter view of all frames
for i = 1:frames
A1(i,:) = -20*log10(abs(fft(ak1(i,:),Ndft)));
A2(i,:) = -20*log10(abs(fft(ak2(i,:),Ndft)));
sd(i) = sum((A1(i,:) - A2(i,:)).^2)/Ndft;
end
printf("sd av %3.2f dB*dB\n", sum(sd)/frames);
figure(1);
clf;
subplot(211)
fs=fopen(raw_filename,"rb");
s = fread(fs,Inf,"short");
plot(s);
subplot(212)
plot(sd);
% now enter single step mode so we can analyse each frame
sn_name = strcat(dump_file_prefix,"_sn.txt");
Sn = load(sn_name);
lsp1_filename = sprintf("%s_lsp.txt", dump_file_prefix);
lsp2_filename = sprintf("%s_lsp_.txt", dump_file_prefix);
lsp1 = load(lsp1_filename);
lsp2 = load(lsp2_filename);
weights_filename = sprintf("%s_weights.txt", dump_file_prefix);
if file_in_path(".",weights_filename)
weights = load(weights_filename);
end
k = ' ';
do
figure(2);
clf;
subplot(211)
s = [ Sn(2*f-1,:) Sn(2*f,:) ];
size(s);
plot(s);
axis([1 length(s) -20000 20000]);
subplot(212);
plot((1:Ndft/2)*4000/(Ndft/2), A1(f,1:(Ndft/2)),";A1;r");
axis([1 4000 -20 40]);
hold on;
plot((1:Ndft/2)*4000/(Ndft/2), A2(f,1:(Ndft/2)),";A2;");
if file_in_path(".",weights_filename)
plot(lsp1(f,:)*4000/pi, weights(f,:),";weights;g+");
end
for l=1:10
plot([lsp1(f,l)*4000/pi lsp1(f,l)*4000/pi], [0 -10], 'r');
plot([lsp2(f,l)*4000/pi lsp2(f,l)*4000/pi], [-10 -20], 'b');
endfor
plot(0,0,';lsp1;r');
plot(0,0,';lsp2;b');
sd_str = sprintf(";sd %3.2f dB*dB;", sd(f));
plot(0,0,sd_str);
hold off;
% interactive menu
printf("\rframe: %d menu: n-next b-back q-quit", f);
fflush(stdout);
k = kbhit();
if (k == 'n')
f = f + 1;
endif
if (k == 'b')
f = f - 1;
endif
until (k == 'q')
printf("\n");
endfunction

86
libs/libcodec2/octave/spec.m Normal file
View File

@ -0,0 +1,86 @@
% spec.m
% Jean Marc Valin
%
% Spectrogram function for Octave
%
% Copyright (c) John-Marc Valin 2012
%
% 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.
%
% - Neither the name of Jean Marc Valin nor the names of its
% contributors may be used to endorse or promote products derived from
% this software without specific prior written permission.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
% ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
% LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
% A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE FOUNDATION OR
% CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
% EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
% PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
% PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
% LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
% NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
% SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
function X = spec(x, Fs, framesize, start, stop)
gr=[zeros(1024,1),[0:1023]'/1023,.68*[0:1023]'/1023];
%gr=[.4*[0:1023]'/1023,[0:1023]'/1023,.68*[0:1023]'/1023];
%t=[0:1023]'/1023;
%t=(1+.25*t-cos(pi*t))/2.25;
%gr = [zeros(1024,1),t,.68*t];
%colormap(gr);
if nargin < 2 || isempty(Fs)
Fs = 44100;
end
if nargin < 3 || isempty(framesize)
framesize = 2048;
endif
offset = framesize/4;
X = 20*log10(abs(specgram(x, framesize, 48000, blackmanharris(framesize)', framesize-offset)));
XM=max(max(X));
X = max(XM-120,X);
%size(X)
F = -[framesize/2-1:-1:0]/framesize*Fs;
%F = [0:127]/128*24000;
T=[1:size(X,2)]/Fs*offset;
%imagesc(X(end:-1:1,:));
if nargin < 4 || isempty(start)
istart=1;
else
istart = round(start*Fs/offset);
end
if nargin < 5 || isempty(stop)
istop = size(X,2);
else
istop = round(stop*Fs/offset);
endif
istart = max(1,istart);
istop = min(istop, size(X,2));
imagesc(T(1+istart:istop), F, X(end:-1:1,1+istart:istop));
X = X(:,1+istart:istop);

306
libs/libcodec2/octave/tfdmdv.m Normal file
View File

@ -0,0 +1,306 @@
% tfdmdv.m
%
% Octave script that tests the C port of the FDMDV modem. This script loads
% the output of unittest/tfdmdv.c and compares it to the output of the
% reference versions of the same functions written in Octave.
%
% Copyright David Rowe 2012
% This program is distributed under the terms of the GNU General Public License
% Version 2
%
fdmdv; % load modem code
% Generate reference vectors using Octave implementation of FDMDV modem
global passes;
global fails;
passes = fails = 0;
frames = 25;
prev_tx_symbols = ones(Nc+1,1);
prev_rx_symbols = ones(Nc+1,1);
foff_phase_rect = 1;
coarse_fine = 0;
fest_state = 0;
channel = [];
channel_count = 0;
next_nin = M;
sig_est = zeros(Nc+1,1);
noise_est = zeros(Nc+1,1);
% Octave outputs we want to collect for comparison to C version
tx_bits_log = [];
tx_symbols_log = [];
tx_baseband_log = [];
tx_fdm_log = [];
pilot_baseband1_log = [];
pilot_baseband2_log = [];
pilot_lpf1_log = [];
pilot_lpf2_log = [];
S1_log = [];
S2_log = [];
foff_coarse_log = [];
foff_fine_log = [];
foff_log = [];
rx_baseband_log = [];
rx_filt_log = [];
env_log = [];
rx_timing_log = [];
rx_symbols_log = [];
rx_bits_log = [];
sync_bit_log = [];
coarse_fine_log = [];
nin_log = [];
sig_est_log = [];
noise_est_log = [];
for f=1:frames
% modulator
tx_bits = get_test_bits(Nc*Nb);
tx_bits_log = [tx_bits_log tx_bits];
tx_symbols = bits_to_qpsk(prev_tx_symbols, tx_bits, 'dqpsk');
prev_tx_symbols = tx_symbols;
tx_symbols_log = [tx_symbols_log tx_symbols];
tx_baseband = tx_filter(tx_symbols);
tx_baseband_log = [tx_baseband_log tx_baseband];
tx_fdm = fdm_upconvert(tx_baseband);
tx_fdm_log = [tx_fdm_log tx_fdm];
% channel
nin = next_nin;
%nin = 120;
%nin = M;
%if (f == 3)
% nin = 120;
%elseif (f == 4)
% nin = 200;
%else
% nin = M;
%end
channel = [channel real(tx_fdm)];
channel_count += M;
rx_fdm = channel(1:nin);
channel = channel(nin+1:channel_count);
channel_count -= nin;
% demodulator
[pilot prev_pilot pilot_lut_index prev_pilot_lut_index] = get_pilot(pilot_lut_index, prev_pilot_lut_index, nin);
[foff_coarse S1 S2] = rx_est_freq_offset(rx_fdm, pilot, prev_pilot, nin);
if coarse_fine == 0
foff = foff_coarse;
end
foff_coarse_log = [foff_coarse_log foff_coarse];
pilot_baseband1_log = [pilot_baseband1_log pilot_baseband1];
pilot_baseband2_log = [pilot_baseband2_log pilot_baseband2];
pilot_lpf1_log = [pilot_lpf1_log pilot_lpf1];
pilot_lpf2_log = [pilot_lpf2_log pilot_lpf2];
S1_log = [S1_log S1];
S2_log = [S2_log S2];
foff_rect = exp(j*2*pi*foff/Fs);
for i=1:nin
foff_phase_rect *= foff_rect';
rx_fdm_fcorr(i) = rx_fdm(i)*foff_phase_rect;
end
rx_baseband = fdm_downconvert(rx_fdm_fcorr, nin);
rx_baseband_log = [rx_baseband_log rx_baseband];
rx_filt = rx_filter(rx_baseband, nin);
rx_filt_log = [rx_filt_log rx_filt];
[rx_symbols rx_timing env] = rx_est_timing(rx_filt, rx_baseband, nin);
env_log = [env_log env];
rx_timing_log = [rx_timing_log rx_timing];
rx_symbols_log = [rx_symbols_log rx_symbols];
next_nin = M;
if rx_timing > 2*M/P
next_nin += M/P;
end
if rx_timing < 0;
next_nin -= M/P;
end
nin_log = [nin_log nin];
[rx_bits sync_bit foff_fine pd] = qpsk_to_bits(prev_rx_symbols, rx_symbols, 'dqpsk');
[sig_est noise_est] = snr_update(sig_est, noise_est, pd);
sig_est_log = [sig_est_log sig_est];
noise_est_log = [noise_est_log noise_est];
prev_rx_symbols = rx_symbols;
rx_bits_log = [rx_bits_log rx_bits];
foff_fine_log = [foff_fine_log foff_fine];
sync_bit_log = [sync_bit_log sync_bit];
foff -= 0.5*foff_fine;
foff_log = [foff_log foff];
% freq est state machine
[coarse_fine fest_state] = freq_state(sync_bit, fest_state);
coarse_fine_log = [coarse_fine_log coarse_fine];
end
% Compare to the output from the C version
load ../unittest/tfdmdv_out.txt
% Helper functions to plot output of C verson and difference between Octave and C versions
function stem_sig_and_error(plotnum, subplotnum, sig, error, titlestr, axisvec)
figure(plotnum)
subplot(subplotnum)
stem(sig);
hold on;
stem(error,'g');
hold off;
if nargin == 6
axis(axisvec);
end
title(titlestr);
endfunction
function plot_sig_and_error(plotnum, subplotnum, sig, error, titlestr, axisvec)
figure(plotnum)
subplot(subplotnum)
plot(sig);
hold on;
plot(error,'g');
hold off;
if nargin == 6
axis(axisvec);
end
title(titlestr);
endfunction
% ---------------------------------------------------------------------------------------
% Plot output and test each C function
% ---------------------------------------------------------------------------------------
% fdmdv_get_test_bits() & bits_to_dqpsk_symbols()
n = 28;
stem_sig_and_error(1, 211, tx_bits_log_c(1:n), tx_bits_log(1:n) - tx_bits_log_c(1:n), 'tx bits', [1 n -1.5 1.5])
stem_sig_and_error(1, 212, real(tx_symbols_log_c(1:n/2)), real(tx_symbols_log(1:n/2) - tx_symbols_log_c(1:n/2)), 'tx symbols real', [1 n/2 -1.5 1.5])
% tx_filter()
diff = tx_baseband_log - tx_baseband_log_c;
c=15;
plot_sig_and_error(2, 211, real(tx_baseband_log_c(c,:)), real(sum(diff)), 'tx baseband real')
plot_sig_and_error(2, 212, imag(tx_baseband_log_c(c,:)), imag(sum(diff)), 'tx baseband imag')
% fdm_upconvert()
plot_sig_and_error(3, 211, real(tx_fdm_log_c), real(tx_fdm_log - tx_fdm_log_c), 'tx fdm real')
plot_sig_and_error(3, 212, imag(tx_fdm_log_c), imag(tx_fdm_log - tx_fdm_log_c), 'tx fdm imag')
% generate_pilot_lut()
plot_sig_and_error(4, 211, real(pilot_lut_c), real(pilot_lut - pilot_lut_c), 'pilot lut real')
plot_sig_and_error(4, 212, imag(pilot_lut_c), imag(pilot_lut - pilot_lut_c), 'pilot lut imag')
% rx_est_freq_offset()
st=1; en = 3*Npilotbaseband;
plot_sig_and_error(5, 211, real(pilot_baseband1_log(st:en)), real(pilot_baseband1_log(st:en) - pilot_baseband1_log_c(st:en)), 'pilot baseband1 real' )
plot_sig_and_error(5, 212, real(pilot_baseband2_log(st:en)), real(pilot_baseband2_log(st:en) - pilot_baseband2_log_c(st:en)), 'pilot baseband2 real' )
st=1; en = 3*Npilotlpf;
plot_sig_and_error(6, 211, real(pilot_lpf1_log(st:en)), real(pilot_lpf1_log(st:en) - pilot_lpf1_log_c(st:en)), 'pilot lpf1 real' )
plot_sig_and_error(6, 212, real(pilot_lpf2_log(st:en)), real(pilot_lpf2_log(st:en) - pilot_lpf2_log_c(st:en)), 'pilot lpf2 real' )
plot_sig_and_error(7, 211, real(S1_log), real(S1_log - S1_log_c), 'S1 real' )
plot_sig_and_error(7, 212, imag(S1_log), imag(S1_log - S1_log_c), 'S1 imag' )
plot_sig_and_error(8, 211, real(S2_log), real(S2_log - S2_log_c), 'S2 real' )
plot_sig_and_error(8, 212, imag(S2_log), imag(S2_log - S2_log_c), 'S2 imag' )
plot_sig_and_error(9, 211, foff_coarse_log, foff_coarse_log - foff_coarse_log_c, 'Coarse Freq Offset' )
plot_sig_and_error(9, 212, foff_fine_log, foff_fine_log - foff_fine_log_c, 'Fine Freq Offset' )
plot_sig_and_error(10, 211, foff_log, foff_log - foff_log_c, 'Freq Offset' )
plot_sig_and_error(10, 212, coarse_fine_log, coarse_fine_log - coarse_fine_log_c, 'Freq Est Coarse(0) Fine(1)', [1 frames -0.5 1.5] )
c=15;
plot_sig_and_error(11, 211, real(rx_baseband_log(c,:)), real(rx_baseband_log(c,:) - rx_baseband_log_c(c,:)), 'Rx baseband real' )
plot_sig_and_error(11, 212, imag(rx_baseband_log(c,:)), imag(rx_baseband_log(c,:) - rx_baseband_log_c(c,:)), 'Rx baseband imag' )
plot_sig_and_error(12, 211, real(rx_filt_log(c,:)), real(rx_filt_log(c,:) - rx_filt_log_c(c,:)), 'Rx filt real' )
plot_sig_and_error(12, 212, imag(rx_filt_log(c,:)), imag(rx_filt_log(c,:) - rx_filt_log_c(c,:)), 'Rx filt imag' )
st=1; en=3*Nt*P;
plot_sig_and_error(13, 211, env_log(st:en), env_log(st:en) - env_log_c(st:en), 'env' )
stem_sig_and_error(13, 212, real(rx_symbols_log(c,:)), real(rx_symbols_log(c,:) - rx_symbols_log_c(c,:)), 'rx symbols' )
st=10*28;
en = 12*28;
plot_sig_and_error(14, 211, rx_timing_log, rx_timing_log - rx_timing_log_c, 'Rx Timing' )
stem_sig_and_error(14, 212, sync_bit_log_c, sync_bit_log - sync_bit_log_c, 'Sync bit', [1 n -1.5 1.5])
stem_sig_and_error(15, 211, rx_bits_log_c(st:en), rx_bits_log(st:en) - rx_bits_log_c(st:en), 'RX bits', [1 en-st -1.5 1.5])
stem_sig_and_error(15, 212, nin_log_c, nin_log - nin_log_c, 'nin')
c = 1;
plot_sig_and_error(16, 211, sig_est_log(c,:), sig_est_log(c,:) - sig_est_log_c(c,:), 'sig est for SNR' )
plot_sig_and_error(16, 212, noise_est_log(c,:), noise_est_log(c,:) - noise_est_log_c(c,:), 'noise est for SNR' )
% ---------------------------------------------------------------------------------------
% AUTOMATED CHECKS ------------------------------------------
% ---------------------------------------------------------------------------------------
function check(a, b, test_name)
global passes;
global fails;
[m n] = size(a);
printf("%s", test_name);
for i=1:(25-length(test_name))
printf(".");
end
printf(": ");
if sum(abs(a - b))/n < 1E-3
printf("OK\n");
passes++;
else
printf("FAIL\n");
fails++;
end
endfunction
check(tx_bits_log, tx_bits_log_c, 'tx_bits');
check(tx_symbols_log, tx_symbols_log_c, 'tx_symbols');
check(tx_baseband_log, tx_baseband_log_c, 'tx_baseband');
check(tx_fdm_log, tx_fdm_log_c, 'tx_fdm');
check(pilot_lut, pilot_lut_c, 'pilot_lut');
check(pilot_baseband1_log, pilot_baseband1_log_c, 'pilot lpf1');
check(pilot_baseband2_log, pilot_baseband2_log_c, 'pilot lpf2');
check(S1_log, S1_log_c, 'S1');
check(S2_log, S2_log_c, 'S2');
check(foff_coarse_log, foff_coarse_log_c, 'foff_coarse');
check(foff_fine_log, foff_fine_log_c, 'foff_fine');
check(foff_log, foff_log_c, 'foff');
check(rx_baseband_log, rx_baseband_log_c, 'rx baseband');
check(rx_filt_log, rx_filt_log_c, 'rx filt');
check(env_log, env_log_c, 'env');
check(rx_timing_log, rx_timing_log_c, 'rx_timing');
check(rx_symbols_log, rx_symbols_log_c, 'rx_symbols');
check(rx_bits_log, rx_bits_log_c, 'rx bits');
check(sync_bit_log, sync_bit_log_c, 'sync bit');
check(coarse_fine_log, coarse_fine_log_c, 'coarse_fine');
check(nin_log, nin_log_c, 'nin');
check(sig_est_log, sig_est_log_c, 'sig_est');
check(noise_est_log, noise_est_log_c, 'noise_est');
printf("\npasses: %d fails: %d\n", passes, fails);

49
libs/libcodec2/octave/tget_spec.m Normal file
View File

@ -0,0 +1,49 @@
% tget-spec.m
%
% Used in conjunction with src/fdmdv_demod to test the
% fdmdv_get_rx_spectrum() function.
%
% codec2-dev/src$ ./fdmdv_demod fdmdv_mod.raw tmp.c2 dump.txt
% octave:3> tget_spec("../src/dump.txt")
%
% Copyright David Rowe 2012
% This program is distributed under the terms of the GNU General Public License
% Version 2
%
function tfft_log(dumpfilename)
load(dumpfilename);
[rows cols] = size(rx_spec_log_c);
Fs = 8000; low_freq = 0; high_freq = 2500;
res = (Fs/2)/cols;
st_bin = low_freq/res + 1;
en_bin = high_freq/res;
xaxis = (st_bin:en_bin)*res;
f_start = 2; f_end = 100;
beta = 0.1;
av = zeros(f_end, en_bin-st_bin+1);
for r=f_start:f_end
x = (1-beta)*av(r-1,:) + beta*rx_spec_log_c(r,st_bin:en_bin);
av(r,:) = x;
end
% spectrogram (waterfall)
figure(1)
clf;
imagesc(av,[-40 0]);
% animated spectrum display
figure(2)
clf;
for r=f_start:f_end
plot(xaxis, av(r,:))
axis([ low_freq high_freq -40 0])
sleep(0.1)
end
endfunction

52
libs/libcodec2/octave/twotone.m Normal file
View File

@ -0,0 +1,52 @@
% twotone.m
% David Rowe Aug 2012
% Used to experiment with combining phase of two tones
function cbphase
Wo = 100.0*pi/4000;
L = floor(pi/Wo);
phi = zeros(1,L);
% two harmonics
a = 20; b = 21;
% set up phases to whatever
phi(a) = -pi;
phi(b) = -pi/2;
% synthesis the two-tone signal
N = 16000;
Nplot = 250;
s = zeros(1,N);
for m=a:b
s_m = cos(m*Wo*(0:(N-1)) + phi(m));
s = s + s_m;
endfor
% now our theory says that this signal should be the same perceptually
phi_(a) = (phi(a) - phi(b))/2;
phi_(b) = (phi(b) - phi(a))/2;
s_ = zeros(1,N);
for m=a:b
s_m = cos(m*Wo*(0:(N-1)) + phi_(m));
s_ = s_ + s_m;
endfor
% plot them and see if envelope has the same phase, but "carriers"
% have different phase
figure(1);
clf;
subplot(211);
plot(s(1:Nplot));
subplot(212);
plot(s_(1:Nplot),'r');
endfunction

76
libs/libcodec2/octave/twotone1.m Normal file
View File

@ -0,0 +1,76 @@
% twotone1.m
% David Rowe Aug 17 2012
%
% Used to experiment with combining phase of two tones. This version
% sets up a complete synthetic speech signal then tries to combine the
% phase of high frequency tones. Lets see if we can do this and keep perceptual
% nature of signal the same.
function twotone1
% hts1a frame 47
Wo = 0.093168;
L = 33;
%A = [69.626907 460.218536 839.677429 2577.498047 972.647888 712.755066 489.048553 364.830536 409.230652 371.767487 489.112854 893.127014 2447.596680 752.878113 475.720520 234.452271 248.161606 232.171051 202.669891 323.914490 678.749451 362.958038 211.652512 170.764435 148.631790 169.261673 272.254150 176.872375 67.344391 99.022301 60.812035 34.319073 14.864757];
A = zeros(1,L)*100;
phi = [1.560274 1.508063 -1.565184 1.289117 -2.547365 1.412528 -1.303992 3.121130 1.087573 -1.158161 -2.928007 0.995093 -2.614023 0.246136 -2.267406 2.143802 -0.273431 -2.266897 1.685171 -0.668712 2.699722 -1.151891 2.406379 -0.046192 -2.718611 0.761067 -2.305014 0.133172 -1.428978 1.492630 -1.668385 1.539734 -1.336615];
%phi = zeros(1,L);
st = floor(L/2);
%st = 1;
A(st:st+5) = 1000;
% now set up phase of signal with phase of upper frequency harmonic
% pairs combined
phi_ = phi;
for m=floor(L/2):2:L
phi_(m) = (phi(m) - phi(m+1))/2;
phi_(m+1) = (phi(m+1) - phi(m))/2;
%phi_(m+1) = 0;
end
% synthesise the signals
N = 16000;
Nplot = 250;
s = zeros(1,N);
for m=st:L
s_m = A(m)*cos(m*Wo*(0:(N-1)) + phi(m));
s = s + s_m;
endfor
s_ = zeros(1,N);
for m=st:L
s_m = A(m)*cos(m*Wo*(0:(N-1)) + phi_(m));
s_ = s_ + s_m;
endfor
% plot them, expect to see similar time domain waveforms
figure(1);
clf;
subplot(211);
plot(s(1:Nplot));
subplot(212);
plot(s_(1:Nplot),'r');
figure(2);
clf;
subplot(211);
plot(s(1:Nplot)-s_(1:Nplot));
% save to disk
gain = 1;
fs=fopen("twotone1_orig.raw","wb");
fwrite(fs,gain*s,"short");
fclose(fs);
fs=fopen("twotone1_comb.raw","wb");
fwrite(fs,gain*s_,"short");
fclose(fs);
endfunction

26
libs/libcodec2/octave/undersample.m Normal file
View File

@ -0,0 +1,26 @@
% undersample.m
% David Rowe 5 Nov 2012
% Testing algorithms for plotting undersampled data for fdmdv2 waveform displays
fs=fopen("../raw/hts1a.raw","rb");
s = fread(fs,Inf,"short");
Fs1=8000;
Fs2=200;
M=Fs1/Fs2;
samples=length(s)/M;
s1=zeros(1,2*samples);
for b=1:samples
st = (b-1)*M + 1;
en = b*M;
s1(2*b-1) = max(s(st:en));
s1(2*b) = min(s(st:en));
end
subplot(211)
plot(s)
subplot(212)
plot(s1);

23
libs/libcodec2/portaudio/Makefile Normal file
View File

@ -0,0 +1,23 @@
# Makefile for Portaudio test programs
# Requires Portaudio V19
CFLAGS = -g -Wall -I../src
LIBS = -lm -lportaudio -pthread
SRC = ../src/fdmdv.c ../src/kiss_fft.c ../src/fifo.c
all: pa_rec pa_play pa_recplay pa_impresp
pa_rec: Makefile pa_rec.c $(SRC)
gcc $(CFLAGS) pa_rec.c $(SRC) -o pa_rec $(LIBS)
pa_play: Makefile pa_play.c $(SRC)
gcc $(CFLAGS) pa_play.c $(SRC) -o pa_play $(LIBS)
pa_recplay: Makefile pa_recplay.c $(SRC)
gcc $(CFLAGS) pa_recplay.c $(SRC) -o pa_recplay $(LIBS)
pa_impresp: Makefile pa_impresp.c $(SRC)
gcc $(CFLAGS) pa_impresp.c $(SRC) -o pa_impresp $(LIBS)
clean:
rm -f pa_rec pa_play pa_recplay

254
libs/libcodec2/portaudio/pa_impresp.c Normal file
View File

@ -0,0 +1,254 @@
/*
pa_impresp.c
David Rowe
August 29 2012
Measures the impulse reponse of the path between the speaker and
microphone. Used to explore why Codec audio quality is
different through a speaker and headphones.
Modified from pa_playrec.c
*/
/*
* $Id: paex_record.c 1752 2011-09-08 03:21:55Z philburk $
*
* This program uses the PortAudio Portable Audio Library.
* For more information see: http://www.portaudio.com
* Copyright (c) 1999-2000 Ross Bencina and Phil Burk
*
* Permission is hereby granted, free of charge, to any person obtaining
* a copy of this software and associated documentation files
* (the "Software"), to deal in the Software without restriction,
* including without limitation the rights to use, copy, modify, merge,
* publish, distribute, sublicense, and/or sell copies of the Software,
* and to permit persons to whom the Software is furnished to do so,
* subject to the following conditions:
*
* The above copyright notice and this permission notice shall be
* included in all copies or substantial portions of the Software.
*
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
* MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
* IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR
* ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF
* CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
* WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
*/
#include <assert.h>
#include <stdio.h>
#include <stdlib.h>
#include "portaudio.h"
#include "fdmdv.h"
#define SAMPLE_RATE 48000 /* 48 kHz sampling rate rec. as we
can trust accuracy of sound
card */
#define N8 160 /* processing buffer size at 8 kHz */
#define N48 (N8*FDMDV_OS) /* processing buffer size at 48 kHz */
#define MEM8 (FDMDV_OS_TAPS/FDMDV_OS)
#define NUM_CHANNELS 2 /* I think most sound cards prefer
stereo, we will convert to mono
as we sample */
#define IMPULSE_AMP 16384 /* amplitide of impulse */
#define IMPULSE_PERIOD 0.1 /* period (dly between impulses) in secs */
/* state information passed to call back */
typedef struct {
float in48k[FDMDV_OS_TAPS + N48];
float in8k[MEM8 + N8];
FILE *fimp;
float *impulse_buf;
int impulse_buf_length;
int impulse_sample_count;
int framesLeft;
} paTestData;
/*
This routine will be called by the PortAudio engine when audio is
required. It may be called at interrupt level on some machines so
don't do anything that could mess up the system like calling
malloc() or free().
*/
static int callback( const void *inputBuffer, void *outputBuffer,
unsigned long framesPerBuffer,
const PaStreamCallbackTimeInfo* timeInfo,
PaStreamCallbackFlags statusFlags,
void *userData )
{
paTestData *data = (paTestData*)userData;
int i;
short *rptr = (short*)inputBuffer;
short *wptr = (short*)outputBuffer;
float *in8k = data->in8k;
float *in48k = data->in48k;
float out8k[N8];
float out48k[N48];
short out48k_short[N48];
short out8k_short[N8];
(void) timeInfo;
(void) statusFlags;
assert(inputBuffer != NULL);
/* just use left channel */
for(i=0; i<framesPerBuffer; i++,rptr+=2)
data->in48k[i+FDMDV_OS_TAPS] = *rptr;
/* downsample and update filter memory */
fdmdv_48_to_8(out8k, &in48k[FDMDV_OS_TAPS], N8);
for(i=0; i<FDMDV_OS_TAPS; i++)
in48k[i] = in48k[i+framesPerBuffer];
/* write impulse response to disk */
for(i=0; i<N8; i++)
out8k_short[i] = out8k[i];
fwrite(out8k_short, sizeof(short), N8, data->fimp);
/* play side, read from impulse buffer */
for(i=0; i<N8; i++) {
in8k[MEM8+i] = data->impulse_buf[data->impulse_sample_count];
data->impulse_sample_count++;
if (data->impulse_sample_count == data->impulse_buf_length)
data->impulse_sample_count = 0;
}
/* upsample and update filter memory */
fdmdv_8_to_48(out48k, &in8k[MEM8], N8);
for(i=0; i<MEM8; i++)
in8k[i] = in8k[i+N8];
assert(outputBuffer != NULL);
/* write signal to both channels */
for(i=0; i<N48; i++)
out48k_short[i] = (short)out48k[i];
for(i=0; i<framesPerBuffer; i++,wptr+=2) {
wptr[0] = out48k_short[i];
wptr[1] = out48k_short[i];
}
data->framesLeft -= framesPerBuffer;
if (data->framesLeft > 0)
return paContinue;
else
return paComplete;
}
int main(int argc, char *argv[])
{
PaStreamParameters inputParameters, outputParameters;
PaStream* stream;
PaError err = paNoError;
paTestData data;
int i, numSecs;
if (argc != 3) {
printf("usage: %s impulseRawFile time(s)\n", argv[0]);
exit(0);
}
data.fimp = fopen(argv[1], "wb");
if (data.fimp == NULL) {
printf("Error opening impulse output file %s\n", argv[1]);
exit(1);
}
numSecs = atoi(argv[2]);
data.framesLeft = numSecs * SAMPLE_RATE;
/* init filter states */
for(i=0; i<MEM8; i++)
data.in8k[i] = 0.0;
for(i=0; i<FDMDV_OS_TAPS; i++)
data.in48k[i] = 0.0;
/* init imupulse */
data.impulse_buf_length = IMPULSE_PERIOD*(SAMPLE_RATE/FDMDV_OS);
printf("%d\n",data.impulse_buf_length);
data.impulse_buf = (float*)malloc(data.impulse_buf_length*sizeof(float));
assert(data.impulse_buf != NULL);
data.impulse_buf[0] = IMPULSE_AMP;
for(i=1; i<data.impulse_buf_length; i++)
data.impulse_buf[i] = 0;
data.impulse_sample_count = 0;
err = Pa_Initialize();
if( err != paNoError ) goto done;
inputParameters.device = Pa_GetDefaultInputDevice(); /* default input device */
if (inputParameters.device == paNoDevice) {
fprintf(stderr,"Error: No default input device.\n");
goto done;
}
inputParameters.channelCount = NUM_CHANNELS; /* stereo input */
inputParameters.sampleFormat = paInt16;
inputParameters.suggestedLatency = Pa_GetDeviceInfo( inputParameters.device )->defaultLowInputLatency;
inputParameters.hostApiSpecificStreamInfo = NULL;
outputParameters.device = Pa_GetDefaultOutputDevice(); /* default output device */
if (outputParameters.device == paNoDevice) {
fprintf(stderr,"Error: No default output device.\n");
goto done;
}
outputParameters.channelCount = NUM_CHANNELS; /* stereo output */
outputParameters.sampleFormat = paInt16;
outputParameters.suggestedLatency = Pa_GetDeviceInfo( outputParameters.device )->defaultLowOutputLatency;
outputParameters.hostApiSpecificStreamInfo = NULL;
/* Play some audio --------------------------------------------- */
err = Pa_OpenStream(
&stream,
&inputParameters,
&outputParameters,
SAMPLE_RATE,
N48,
paClipOff,
callback,
&data );
if( err != paNoError ) goto done;
err = Pa_StartStream( stream );
if( err != paNoError ) goto done;
while( ( err = Pa_IsStreamActive( stream ) ) == 1 )
{
Pa_Sleep(100);
}
if( err < 0 ) goto done;
err = Pa_CloseStream( stream );
if( err != paNoError ) goto done;
done:
Pa_Terminate();
if( err != paNoError )
{
fprintf( stderr, "An error occured while using the portaudio stream\n" );
fprintf( stderr, "Error number: %d\n", err );
fprintf( stderr, "Error message: %s\n", Pa_GetErrorText( err ) );
err = 1; /* Always return 0 or 1, but no other return codes. */
}
fclose(data.fimp);
return err;
}

199
libs/libcodec2/portaudio/pa_play.c Normal file
View File

@ -0,0 +1,199 @@
/*
pa_play.c
David Rowe
July 8 2012
Converts samples from a 16 bit short 8000 Hz rawfile to 480000Hz
sample rate and plays them using the default sound device. Used as
an intermediate step in Portaudio integration.
Modified from paex_record.c Portaudio example. Original author
author Phil Burk http://www.softsynth.com
*/
/*
* $Id: paex_record.c 1752 2011-09-08 03:21:55Z philburk $
*
* This program uses the PortAudio Portable Audio Library.
* For more information see: http://www.portaudio.com
* Copyright (c) 1999-2000 Ross Bencina and Phil Burk
*
* Permission is hereby granted, free of charge, to any person obtaining
* a copy of this software and associated documentation files
* (the "Software"), to deal in the Software without restriction,
* including without limitation the rights to use, copy, modify, merge,
* publish, distribute, sublicense, and/or sell copies of the Software,
* and to permit persons to whom the Software is furnished to do so,
* subject to the following conditions:
*
* The above copyright notice and this permission notice shall be
* included in all copies or substantial portions of the Software.
*
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
* MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
* IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR
* ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF
* CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
* WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
*/
#include <assert.h>
#include <stdio.h>
#include <stdlib.h>
#include "portaudio.h"
#include "fdmdv.h"
#define SAMPLE_RATE 48000 /* 48 kHz sampling rate rec. as we
can trust accuracy of sound
card */
#define N8 160 /* processing buffer size at 8 kHz */
#define N48 (N8*FDMDV_OS) /* processing buffer size at 48 kHz */
#define MEM8 (FDMDV_OS_TAPS/FDMDV_OS)
#define NUM_CHANNELS 2 /* I think most sound cards prefer
stereo, we will convert to mono
as we sample */
/* state information passed to call back */
typedef struct {
FILE *fin;
float in8k[MEM8 + N8];
} paTestData;
/*
This routine will be called by the PortAudio engine when audio is
required. It may be called at interrupt level on some machines so
don't do anything that could mess up the system like calling
malloc() or free().
*/
static int playCallback( const void *inputBuffer, void *outputBuffer,
unsigned long framesPerBuffer,
const PaStreamCallbackTimeInfo* timeInfo,
PaStreamCallbackFlags statusFlags,
void *userData )
{
paTestData *data = (paTestData*)userData;
FILE *fin = data->fin;
int i, nread;
int finished;
short *wptr = (short*)outputBuffer;
float *in8k = data->in8k;
float out48k[N48];
short out48k_short[N48];
short in8k_short[N8];
(void) outputBuffer; /* Prevent unused variable warnings. */
(void) timeInfo;
(void) statusFlags;
(void) userData;
/* note Portaudio docs recs. against making systems calls like
fwrite() in this callback but seems to work OK */
nread = fread(in8k_short, sizeof(short), N8, fin);
if (nread == N8)
finished = paContinue;
else
finished = paComplete;
for(i=0; i<N8; i++)
in8k[MEM8+i] = in8k_short[i];
/* upsample and update filter memory */
fdmdv_8_to_48(out48k, &in8k[MEM8], N8);
for(i=0; i<MEM8; i++)
in8k[i] = in8k[i+N8];
assert(outputBuffer != NULL);
/* write signal to both channels */
for(i=0; i<N48; i++)
out48k_short[i] = (short)out48k[i];
for(i=0; i<framesPerBuffer; i++,wptr+=2) {
wptr[0] = out48k_short[i];
wptr[1] = out48k_short[i];
}
return finished;
}
int main(int argc, char *argv[])
{
PaStreamParameters outputParameters;
PaStream* stream;
PaError err = paNoError;
paTestData data;
int i;
if (argc != 2) {
printf("usage: %s rawFile\n", argv[0]);
exit(0);
}
data.fin = fopen(argv[1], "rt");
if (data.fin == NULL) {
printf("Error opening input raw file %s\n", argv[1]);
exit(1);
}
for(i=0; i<MEM8; i++)
data.in8k[i] = 0.0;
err = Pa_Initialize();
if( err != paNoError ) goto done;
outputParameters.device = Pa_GetDefaultOutputDevice(); /* default input device */
if (outputParameters.device == paNoDevice) {
fprintf(stderr,"Error: No default output device.\n");
goto done;
}
outputParameters.channelCount = NUM_CHANNELS; /* stereo input */
outputParameters.sampleFormat = paInt16;
outputParameters.suggestedLatency = Pa_GetDeviceInfo( outputParameters.device )->defaultLowOutputLatency;
outputParameters.hostApiSpecificStreamInfo = NULL;
/* Play some audio --------------------------------------------- */
err = Pa_OpenStream(
&stream,
NULL,
&outputParameters,
SAMPLE_RATE,
N48,
paClipOff,
playCallback,
&data );
if( err != paNoError ) goto done;
err = Pa_StartStream( stream );
if( err != paNoError ) goto done;
while( ( err = Pa_IsStreamActive( stream ) ) == 1 )
{
Pa_Sleep(100);
}
if( err < 0 ) goto done;
err = Pa_CloseStream( stream );
if( err != paNoError ) goto done;
fclose(data.fin);
done:
Pa_Terminate();
if( err != paNoError )
{
fprintf( stderr, "An error occured while using the portaudio stream\n" );
fprintf( stderr, "Error number: %d\n", err );
fprintf( stderr, "Error message: %s\n", Pa_GetErrorText( err ) );
err = 1; /* Always return 0 or 1, but no other return codes. */
}
return err;
}

212
libs/libcodec2/portaudio/pa_rec.c Normal file
View File

@ -0,0 +1,212 @@
/*
pa_rec.c
David Rowe
July 6 2012
Records at 48000 Hz from default sound device, convertes to 8 kHz,
and saves to raw file. Used to get experience with Portaudio.
Modified from paex_record.c Portaudio example. Original author
author Phil Burk http://www.softsynth.com
To Build:
gcc paex_rec.c -o paex_rec -lm -lrt -lportaudio -pthread
*/
/*
* $Id: paex_record.c 1752 2011-09-08 03:21:55Z philburk $
*
* This program uses the PortAudio Portable Audio Library.
* For more information see: http://www.portaudio.com
* Copyright (c) 1999-2000 Ross Bencina and Phil Burk
*
* Permission is hereby granted, free of charge, to any person obtaining
* a copy of this software and associated documentation files
* (the "Software"), to deal in the Software without restriction,
* including without limitation the rights to use, copy, modify, merge,
* publish, distribute, sublicense, and/or sell copies of the Software,
* and to permit persons to whom the Software is furnished to do so,
* subject to the following conditions:
*
* The above copyright notice and this permission notice shall be
* included in all copies or substantial portions of the Software.
*
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
* MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
* IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR
* ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF
* CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
* WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
*/
#include <assert.h>
#include <stdio.h>
#include <stdlib.h>
#include "portaudio.h"
#include "fdmdv.h"
#define SAMPLE_RATE 48000 /* 48 kHz sampling rate rec. as we
can trust accuracy of sound
card */
#define N8 160 /* processing buffer size at 8 kHz */
#define N48 (N8*FDMDV_OS) /* processing buffer size at 48 kHz */
#define NUM_CHANNELS 2 /* I think most sound cards prefer
stereo, we will convert to mono
as we sample */
/* state information passed to call back */
typedef struct {
FILE *fout;
int framesLeft;
float in48k[FDMDV_OS_TAPS + N48];
} paTestData;
/*
This routine will be called by the PortAudio engine when audio is
available. It may be called at interrupt level on some machines so
don't do anything that could mess up the system like calling
malloc() or free().
*/
static int recordCallback( const void *inputBuffer, void *outputBuffer,
unsigned long framesPerBuffer,
const PaStreamCallbackTimeInfo* timeInfo,
PaStreamCallbackFlags statusFlags,
void *userData )
{
paTestData *data = (paTestData*)userData;
FILE *fout = data->fout;
int framesToCopy;
int i;
int finished;
short *rptr = (short*)inputBuffer;
float out8k[N8];
short out8k_short[N8];
(void) outputBuffer; /* Prevent unused variable warnings. */
(void) timeInfo;
(void) statusFlags;
(void) userData;
if (data->framesLeft < framesPerBuffer) {
framesToCopy = data->framesLeft;
finished = paComplete;
}
else {
framesToCopy = framesPerBuffer;
finished = paContinue;
}
data->framesLeft -= framesToCopy;
assert(inputBuffer != NULL);
/* just use left channel */
for(i=0; i<framesToCopy; i++,rptr+=2)
data->in48k[i+FDMDV_OS_TAPS] = *rptr;
/* downsample and update filter memory */
fdmdv_48_to_8(out8k, &data->in48k[FDMDV_OS_TAPS], N8);
for(i=0; i<FDMDV_OS_TAPS; i++)
data->in48k[i] = data->in48k[i+framesToCopy];
/* save 8k to disk */
for(i=0; i<N8; i++)
out8k_short[i] = (short)out8k[i];
/* note Portaudio docs recs. against making systems calls like
fwrite() in this callback but seems to work OK */
fwrite(out8k_short, sizeof(short), N8, fout);
return finished;
}
int main(int argc, char *argv[])
{
PaStreamParameters inputParameters;
PaStream* stream;
PaError err = paNoError;
paTestData data;
int i;
int numSecs;
if (argc != 3) {
printf("usage: %s rawFile time(s)\n", argv[0]);
exit(0);
}
data.fout = fopen(argv[1], "wt");
if (data.fout == NULL) {
printf("Error opening output raw file %s\n", argv[1]);
exit(1);
}
numSecs = atoi(argv[2]);
data.framesLeft = numSecs * SAMPLE_RATE;
for(i=0; i<FDMDV_OS_TAPS; i++)
data.in48k[i] = 0.0;
err = Pa_Initialize();
if( err != paNoError ) goto done;
printf( "PortAudio version number = %d\nPortAudio version text = '%s'\n",
Pa_GetVersion(), Pa_GetVersionText() );
inputParameters.device = Pa_GetDefaultInputDevice(); /* default input device */
if (inputParameters.device == paNoDevice) {
fprintf(stderr,"Error: No default input device.\n");
goto done;
}
inputParameters.channelCount = NUM_CHANNELS; /* stereo input */
inputParameters.sampleFormat = paInt16;
inputParameters.suggestedLatency = Pa_GetDeviceInfo( inputParameters.device )->defaultLowInputLatency;
inputParameters.hostApiSpecificStreamInfo = NULL;
/* Record some audio --------------------------------------------- */
err = Pa_OpenStream(
&stream,
&inputParameters,
NULL, /* &outputParameters, */
SAMPLE_RATE,
N48,
paClipOff, /* we won't output out of range samples so don't bother clipping them */
recordCallback,
&data );
if( err != paNoError ) goto done;
err = Pa_StartStream( stream );
if( err != paNoError ) goto done;
while( ( err = Pa_IsStreamActive( stream ) ) == 1 )
{
Pa_Sleep(100);
}
if( err < 0 ) goto done;
err = Pa_CloseStream( stream );
if( err != paNoError ) goto done;
fclose(data.fout);
done:
Pa_Terminate();
if( err != paNoError )
{
fprintf( stderr, "An error occured while using the portaudio stream\n" );
fprintf( stderr, "Error number: %d\n", err );
fprintf( stderr, "Error message: %s\n", Pa_GetErrorText( err ) );
err = 1; /* Always return 0 or 1, but no other return codes. */
}
return err;
}

270
libs/libcodec2/portaudio/pa_recplay.c Normal file
View File

@ -0,0 +1,270 @@
/*
pa_recplay.c
David Rowe
July 8 2012
Echos audio from sound card input to sound card output. Samples at
48 kHz, converts to 8 kHz, converts back to 48kHz, and plays using
the default sound device. Used as an intermediate step in
Portaudio integration.
Modified from paex_record.c Portaudio example. Original author
author Phil Burk http://www.softsynth.com
*/
/*
* $Id: paex_record.c 1752 2011-09-08 03:21:55Z philburk $
*
* This program uses the PortAudio Portable Audio Library.
* For more information see: http://www.portaudio.com
* Copyright (c) 1999-2000 Ross Bencina and Phil Burk
*
* Permission is hereby granted, free of charge, to any person obtaining
* a copy of this software and associated documentation files
* (the "Software"), to deal in the Software without restriction,
* including without limitation the rights to use, copy, modify, merge,
* publish, distribute, sublicense, and/or sell copies of the Software,
* and to permit persons to whom the Software is furnished to do so,
* subject to the following conditions:
*
* The above copyright notice and this permission notice shall be
* included in all copies or substantial portions of the Software.
*
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
* MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
* IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR
* ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF
* CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
* WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
*/
#include <assert.h>
#include <stdio.h>
#include <stdlib.h>
#include "portaudio.h"
#include "fdmdv.h"
#include "fifo.h"
#define SAMPLE_RATE 48000 /* 48 kHz sampling rate rec. as we
can trust accuracy of sound
card */
#define N8 160 /* processing buffer size at 8 kHz */
#define N48 (N8*FDMDV_OS) /* processing buffer size at 48 kHz */
#define MEM8 (FDMDV_OS_TAPS/FDMDV_OS)
#define NUM_CHANNELS 2 /* I think most sound cards prefer
stereo, we will convert to mono
as we sample */
#define MAX_FPB 2048 /* maximum value of framesPerBuffer */
/* state information passed to call back */
typedef struct {
float in48k[FDMDV_OS_TAPS + N48];
float in8k[MEM8 + N8];
struct FIFO *infifo;
struct FIFO *outfifo;
} paTestData;
/*
This routine will be called by the PortAudio engine when audio is
required. It may be called at interrupt level on some machines so
don't do anything that could mess up the system like calling
malloc() or free().
*/
static int callback( const void *inputBuffer, void *outputBuffer,
unsigned long framesPerBuffer,
const PaStreamCallbackTimeInfo* timeInfo,
PaStreamCallbackFlags statusFlags,
void *userData )
{
paTestData *data = (paTestData*)userData;
int i;
short *rptr = (short*)inputBuffer;
short *wptr = (short*)outputBuffer;
float *in8k = data->in8k;
float *in48k = data->in48k;
float out8k[N8];
float out48k[N48];
short out48k_short[N48];
short in48k_short[N48];
short indata[MAX_FPB];
short outdata[MAX_FPB];
(void) timeInfo;
(void) statusFlags;
assert(inputBuffer != NULL);
assert(outputBuffer != NULL);
/*
framesPerBuffer is portaudio-speak for number of samples we
actually get from the record side and need to provide to the
play side. On Linux (at least) it was found that
framesPerBuffer may not always be what we ask for in the
framesPerBuffer field of Pa_OpenStream. For example a request
for 960 sample buffers lead to framesPerBuffer = 1024.
To perform the 48 to 8 kHz conversion we need an integer
multiple of FDMDV_OS samples to support the interpolation and
decimation. As we can't guarantee the size of framesPerBuffer
we do a little FIFO buffering.
*/
//printf("framesPerBuffer: %d N48 %d\n", framesPerBuffer, N48);
/* assemble a mono buffer (just use left channel) and write to FIFO */
assert(framesPerBuffer < MAX_FPB);
for(i=0; i<framesPerBuffer; i++,rptr+=2)
indata[i] = *rptr;
fifo_write(data->infifo, indata, framesPerBuffer);
/* while we have enough samples available ... */
//printf("infifo before: %d\n", fifo_n(data->infifo));
while (fifo_read(data->infifo, in48k_short, N48) == 0) {
/* convert to float */
for(i=0; i<N48; i++)
in48k[FDMDV_OS_TAPS + i] = in48k_short[i];
/* downsample and update filter memory */
fdmdv_48_to_8(out8k, &in48k[FDMDV_OS_TAPS], N8);
for(i=0; i<FDMDV_OS_TAPS; i++)
in48k[i] = in48k[i+N48];
/* play side, back up to 8k */
for(i=0; i<N8; i++)
in8k[MEM8+i] = out8k[i];
/* upsample and update filter memory */
fdmdv_8_to_48(out48k, &in8k[MEM8], N8);
for(i=0; i<MEM8; i++)
in8k[i] = in8k[i+N8];
/* write signal to both channels */
for(i=0; i<N48; i++)
out48k_short[i] = (short)out48k[i];
fifo_write(data->outfifo, out48k_short, N48);
}
//printf("infifo after: %d\n", fifo_n(data->infifo));
//printf("outfifo : %d\n", fifo_n(data->outfifo));
/* OK now set up output samples */
if (fifo_read(data->outfifo, outdata, framesPerBuffer) == 0) {
/* write signal to both channels */
for(i=0; i<framesPerBuffer; i++,wptr+=2) {
wptr[0] = outdata[i];
wptr[1] = outdata[i];
}
}
else {
//printf("no data\n");
/* zero output if no data available */
for(i=0; i<framesPerBuffer; i++,wptr+=2) {
wptr[0] = 0;
wptr[1] = 0;
}
}
return paContinue;
}
int main(int argc, char *argv[])
{
PaStreamParameters inputParameters, outputParameters;
PaStream* stream;
PaError err = paNoError;
paTestData data;
int i;
/* init callback data */
for(i=0; i<MEM8; i++)
data.in8k[i] = 0.0;
for(i=0; i<FDMDV_OS_TAPS; i++)
data.in48k[i] = 0.0;
data.infifo = fifo_create(2*N48);
data.outfifo = fifo_create(2*N48);
/* init port audio */
err = Pa_Initialize();
if( err != paNoError ) goto done;
inputParameters.device = Pa_GetDefaultInputDevice(); /* default input device */
if (inputParameters.device == paNoDevice) {
fprintf(stderr,"Error: No default input device.\n");
goto done;
}
inputParameters.channelCount = NUM_CHANNELS; /* stereo input */
inputParameters.sampleFormat = paInt16;
inputParameters.suggestedLatency = Pa_GetDeviceInfo( inputParameters.device )->defaultLowInputLatency;
inputParameters.hostApiSpecificStreamInfo = NULL;
outputParameters.device = Pa_GetDefaultOutputDevice(); /* default output device */
if (outputParameters.device == paNoDevice) {
fprintf(stderr,"Error: No default output device.\n");
goto done;
}
outputParameters.channelCount = NUM_CHANNELS; /* stereo output */
outputParameters.sampleFormat = paInt16;
outputParameters.suggestedLatency = Pa_GetDeviceInfo( outputParameters.device )->defaultLowOutputLatency;
outputParameters.hostApiSpecificStreamInfo = NULL;
/* Play some audio --------------------------------------------- */
err = Pa_OpenStream(
&stream,
&inputParameters,
&outputParameters,
SAMPLE_RATE,
512,
paClipOff,
callback,
&data );
if( err != paNoError ) goto done;
err = Pa_StartStream( stream );
if( err != paNoError ) goto done;
while( ( err = Pa_IsStreamActive( stream ) ) == 1 )
{
Pa_Sleep(100);
}
if( err < 0 ) goto done;
err = Pa_CloseStream( stream );
if( err != paNoError ) goto done;
done:
Pa_Terminate();
if( err != paNoError )
{
fprintf( stderr, "An error occured while using the portaudio stream\n" );
fprintf( stderr, "Error number: %d\n", err );
fprintf( stderr, "Error message: %s\n", Pa_GetErrorText( err ) );
err = 1; /* Always return 0 or 1, but no other return codes. */
}
fifo_destroy(data.infifo);
fifo_destroy(data.outfifo);
return err;
}

BIN
libs/libcodec2/raw/cq_ref.raw Normal file
View File

Binary file not shown.

BIN
libs/libcodec2/raw/cross.raw Normal file
View File

Binary file not shown.

BIN
libs/libcodec2/raw/cross_melp2400.raw Normal file
View File

Binary file not shown.

BIN
libs/libcodec2/raw/forig_ambe2000.raw Normal file
View File

Binary file not shown.

BIN
libs/libcodec2/raw/g3plx.raw Normal file
View File

Binary file not shown.

BIN
libs/libcodec2/raw/hts1a_ambe2000.raw Normal file
View File

Binary file not shown.

BIN
libs/libcodec2/raw/hts1a_melp.raw Normal file
View File

Binary file not shown.

BIN
libs/libcodec2/raw/hts2a_ambe2000.raw Normal file
View File

Binary file not shown.

BIN
libs/libcodec2/raw/hts2a_melp.raw Normal file
View File

Binary file not shown.

BIN
libs/libcodec2/raw/k6hx.raw Normal file
View File

Binary file not shown.

BIN
libs/libcodec2/raw/kristoff.raw Normal file
View File

Binary file not shown.

BIN
libs/libcodec2/raw/mmt1_ambe2000.raw Normal file
View File

Binary file not shown.

BIN
libs/libcodec2/raw/morig_ambe2000.raw Normal file
View File

Binary file not shown.

BIN
libs/libcodec2/src/c2demo Executable file
View File

Binary file not shown.

BIN
libs/libcodec2/src/fdmdv_demod Executable file
View File

Binary file not shown.

BIN
libs/libcodec2/src/fdmdv_get_test_bits Executable file
View File

Binary file not shown.

BIN
libs/libcodec2/src/fdmdv_interleave Executable file
View File

Binary file not shown.

BIN
libs/libcodec2/src/fdmdv_mod Executable file
View File

Binary file not shown.

BIN
libs/libcodec2/src/fdmdv_put_test_bits Executable file
View File

Binary file not shown.

BIN
libs/libcodec2/src/generate_codebook Executable file
View File

Binary file not shown.

BIN
libs/libcodec2/src/genlspdtcb Executable file
View File

Binary file not shown.

43
libs/libcodec2/unittest/README Normal file
View File

@ -0,0 +1,43 @@
README
for codec2/unittest
Created David Rowe 31 July 2012
Training (experimental) sparse phase VQs:
1/ In ../src/phase.c phase_experiment() enable:
print_sparse_pred_error()
and 'make' c2sim (in src)
2/ Run over a training database:
$ ./c2sim /xhome1/codec2/samples/train.spc --phaseexp > train_phtrain.txt
a) check stats in Octave:
octave> load ../src/train_phtrain.txt
octave> std(nonzeros(train_phtrain(:,1:20)))
octave> hist(nonzeros(train_phtrain(:,1:20)),20)
3/ Extract and convert to floats vector you wish to train for example
first 20 (out of MAX_AMP == 80):
$ ./extract ../src/train_phtrain.txt train_phtrain.flt 1 20
4/ Convert to rectangular:
$ ./polar2rect train_phtrain.flt train_phtrainr.flt
5/ Run this program:
$ ./vqtrainph train_phtrainr.flt 20 1024 vq.txt
Ouput is vq.txt
Tests
-----
+ build up insmallest possible stesp
+ impl errors v alg errors
+ use actual phase data as codebook
+ test vq with rand phases first or known data

BIN
libs/libcodec2/unittest/create_interleaver Executable file
View File

Binary file not shown.

48
libs/libcodec2/unittest/create_interleaver.c Normal file
View File

@ -0,0 +1,48 @@
/*
create_interleaver.c
David Rowe
May 27 2012
Creates an interleaver for Codec 2.
*/
#include <assert.h>
#include <stdlib.h>
#include <stdio.h>
int main(int argc, char * argv[]) {
int m,i, src_bit, dest_bit;
FILE *f;
int *interleaver;
if (argc != 3) {
printf("usage: %s InterleaverBits InterleaverFile\n", argv[0]);
exit(1);
}
m = atoi(argv[1]);
f = fopen(argv[2],"wt");
assert(f != NULL);
interleaver = (int*)malloc(m*sizeof(int));
assert(interleaver != NULL);
for(i=0; i<m; i++)
interleaver[i] = -1;
src_bit = 0;
while(src_bit != m) {
dest_bit = ((float)rand()/RAND_MAX)*m;
if (interleaver[dest_bit] == -1) {
interleaver[dest_bit] = src_bit;
src_bit++;
}
}
for(i=0; i<m; i++) {
fprintf(f, "%d\n", interleaver[i]);
}
fclose(f);
return 0;
}

BIN
libs/libcodec2/unittest/de Executable file
View File

Binary file not shown.

59
libs/libcodec2/unittest/de.c Normal file
View File

@ -0,0 +1,59 @@
/*
de.c
David Rowe
Sep 26 2012
Takes audio from a file, de-emphasises, and sends to output file.
*/
#include <assert.h>
#include <errno.h>
#include <math.h>
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include "lpc.h"
#define N 80
int main(int argc, char *argv[]) {
FILE *fin, *fout;
short buf[N];
float Sn[N], Sn_de[N];
float de_mem = 0.0;
int i;
if (argc != 3) {
printf("usage: de InputRawSpeechFile OutputRawSpeechFile\n");
printf("e.g de input.raw output.raw");
exit(1);
}
if (strcmp(argv[1], "-") == 0) fin = stdin;
else if ( (fin = fopen(argv[1],"rb")) == NULL ) {
fprintf(stderr, "Error opening input speech file: %s: %s.\n",
argv[1], strerror(errno));
exit(1);
}
if (strcmp(argv[2], "-") == 0) fout = stdout;
else if ( (fout = fopen(argv[2],"wb")) == NULL ) {
fprintf(stderr, "Error opening output speech file: %s: %s.\n",
argv[2], strerror(errno));
exit(1);
}
while(fread(buf, sizeof(short), N, fin) == N) {
for(i=0; i<N; i++)
Sn[i] = buf[i];
de_emp(Sn_de, Sn, &de_mem, N);
for(i=0; i<N; i++)
buf[i] = Sn_de[i];
fwrite(buf, sizeof(short), N, fout);
}
fclose(fin);
fclose(fout);
return 0;
}

385
libs/libcodec2/unittest/dvdongle2.c Normal file
View File

@ -0,0 +1,385 @@
/*---------------------------------------------------------------------------*\
FILE........: dvdongle2.c
AUTHOR......: David Rowe
DATE CREATED: 28 Oct 2010
Program to encode and decode raw speech samples using the AMBE codec
implemented on a DV Dongle.
The DV Dongle connects to a USB port and provides encoding and
decoding of compressed audio using the DVSI AMBE2000 full duplex
vocoder DSP chip.
Refs:
[1] http://www.dvdongle.com/
[2] http://www.moetronix.com/files/dvdongletechref100.pdf
[3] http://www.dvsinc.com/manuals/AMBE-2000_manual.pdf
[4] http://www.moetronix.com/files/ambetest103.zip
Serial code based on ser.c sample from http://www.captain.at
Compile with:
gcc dvdongle2.c -o dvdongle2 -Wall -g -O2
Note: This program is not very stable, it sometimes stops part way
through processing an utterance. I made it just good enough to work
most of the time, as my purpose was just to process a few sample
files.
\*---------------------------------------------------------------------------*/
/*
Copyright (C) 1990-2010 David Rowe
All rights reserved.
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License version 2.1, as
published by the Free Software Foundation. This program is
distributed in the hope that it will be useful, but WITHOUT ANY
WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
License for more details.
You should have received a copy of the GNU Lesser General Public License
along with this program; if not, see <http://www.gnu.org/licenses/>.
*/
#include <assert.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <unistd.h>
#include <fcntl.h>
#include <errno.h>
#include <termios.h>
#define MAX_STR 1024
#define LEN_TARGET_NAME_RESPONSE 14
#define N 160
/* message parsing state machine states */
#define MSGSTATE_HDR1 0
#define MSGSTATE_HDR2 1
#define MSGSTATE_DATA 2
#define LENGTH_MASK 0x1FFF /* mask for message length */
#define TYPE_MASK 0xE0 /* mask for upper byte of header */
#define TYPE_C 0x20 /* compressed speech from target */
#define TYPE_UC 0x40 /* uncompressed speech from target */
#define MAX_MSG_LEN 8192
/* Control items sent to DV Dongle */
char target_name[] = {0x04, 0x20, 0x01, 0x00};
/* note [2] appears to be in error, specifies run as 0x02, stop as 0x01 */
char run_state_stop[] = {0x05, 0x00, 0x18, 0x00, 0x00};
char run_state_run[] = {0x05, 0x00, 0x18, 0x00, 0x01};
/* Control item codes from DV Dongle */
char data_item_0[] = {0x42, 0x81};
char data_item_1[] = {0x32, 0xa0};
char run_state[] = {0x05, 0x00};
char idle[] = {0x00, 0x00};
typedef struct {
short header;
char power;
char control1;
short rate[5];
short unused[3];
short dtmf;
short control2;
short channel_data[12];
} COMPRESSED;
COMPRESSED c_in;
COMPRESSED c_out;
FILE *fin, *fout, *f;
int fd, c_msg, uc_msg;
int initport(int fd) {
struct termios options;
// Set the options for the port...
cfmakeraw(&options);
cfsetispeed(&options, B230400);
cfsetospeed(&options, B230400);
options.c_cflag |= (CLOCAL | CREAD);
tcsetattr(fd, TCSANOW, &options);
return 1;
}
int getbaud(int fd) {
struct termios termAttr;
int inputSpeed = -1;
speed_t baudRate;
tcgetattr(fd, &termAttr);
/* Get the input speed */
baudRate = cfgetispeed(&termAttr);
switch (baudRate) {
case B0: inputSpeed = 0; break;
case B50: inputSpeed = 50; break;
case B110: inputSpeed = 110; break;
case B134: inputSpeed = 134; break;
case B150: inputSpeed = 150; break;
case B200: inputSpeed = 200; break;
case B300: inputSpeed = 300; break;
case B600: inputSpeed = 600; break;
case B1200: inputSpeed = 1200; break;
case B1800: inputSpeed = 1800; break;
case B2400: inputSpeed = 2400; break;
case B4800: inputSpeed = 4800; break;
case B9600: inputSpeed = 9600; break;
case B19200: inputSpeed = 19200; break;
case B38400: inputSpeed = 38400; break;
case B57600: inputSpeed = 38400; break;
case B115200: inputSpeed = 38400; break;
case B230400: inputSpeed = 230400; break;
}
return inputSpeed;
}
void write_dongle(int fd, char *data, int len) {
int n;
//printf(" writing %d bytes\n", len);
n = write(fd, data, len);
if (n < 0) {
perror("write failed");
exit(1);
}
}
void read_dongle(int fd, char *data, int len) {
int n;
//printf(" reading %d bytes \n", len);
n = read(fd, data, len);
if (n < 0) {
perror("read failed");
exit(1);
}
//printf(" read %d bytes\n", len);
}
void parse_message(int msg_type, int msg_len, char msg_data[]) {
short buf[N];
COMPRESSED *c_out;
//printf("msg_type: 0x%02x msg_len: %d\n", msg_type, msg_len);
/* echo compressed speech frames back to target */
if (msg_type == TYPE_C) {
c_out = (COMPRESSED*)msg_data;
#ifdef TMP
printf("control1 0x%04x\n", c_out->control1 & 0xff);
printf("rate[0] 0x%04x\n", c_out->rate[0]);
printf("rate[1] 0x%04x\n", c_out->rate[1]);
printf("rate[2] 0x%04x\n", c_out->rate[2]);
printf("rate[3] 0x%04x\n", c_out->rate[3]);
printf("rate[4] 0x%04x\n", c_out->rate[4]);
printf("control2 0x%04x\n", c_out->control2 & 0xffff);
printf("cd[0] 0x%04x\n", c_out->channel_data[0] & 0xffff);
printf("cd[1] 0x%04x\n", c_out->channel_data[1] & 0xffff);
printf("cd[2] 0x%04x\n", c_out->channel_data[2] & 0xffff);
printf("cd[3] 0x%04x\n", c_out->channel_data[3] & 0xffff);
printf("cd[4] 0x%04x\n", c_out->channel_data[4] & 0xffff);
printf("cd[5] 0x%04x\n", c_out->channel_data[5] & 0xffff);
printf("cd[6] 0x%04x\n", c_out->channel_data[6] & 0xffff);
printf("uc_msg %d\n", uc_msg);
#endif
printf("bit errors %d\n", c_out->unused[2]);
memcpy(&c_in.channel_data,
&c_out->channel_data,
sizeof(c_in.channel_data));
write_dongle(fd, data_item_1, sizeof(data_item_1));
write_dongle(fd, (char*)&c_in, sizeof(c_in));
c_msg++;
}
/* write speech buffers to disk */
if (msg_type == TYPE_UC) {
if (fout != NULL) {
fwrite(msg_data, sizeof(char), msg_len-2, fout);
printf("msg_len %d\n", msg_len);
}
if (fin != NULL)
fread(buf, sizeof(short), N, fin);
else
memset(buf, 0, sizeof(buf));
write_dongle(fd, data_item_0, sizeof(data_item_0));
write_dongle(fd, (char*)buf, sizeof(buf));
uc_msg++;
}
}
int main(int argc, char **argv) {
char response[MAX_STR];
int i;
int state, next_state;
short header;
int msg_type, msg_length;
char msg_data[MAX_MSG_LEN];
int n, length;
int r;
char data;
f = fopen("/tmp/log.txt", "wt");
assert(f != NULL);
/* open and configure serial port */
fd = open("/dev/ttyUSB0", O_RDWR | O_NOCTTY | O_NDELAY);
if (fd == -1) {
perror("open_port: Unable to open /dev/ttyS0 - ");
exit(1);
} else {
fcntl(fd, F_SETFL, 0);
}
initport(fd);
fin = NULL;
if (argc >= 2) {
fin = fopen(argv[1],"rb");
assert(fin != NULL);
}
fout = NULL;
if (argc == 3) {
fout = fopen(argv[2],"wb");
assert(fout != NULL);
}
/* check DV Dongle is alive */
write_dongle(fd, target_name, sizeof(target_name));
read_dongle(fd, response, LEN_TARGET_NAME_RESPONSE);
if (strcmp(&response[4],"DV Dongle") != 0) {
printf("DV Dongle not responding\n");
exit(1);
}
printf("Found DV Dongle....\n");
c_in.header = 0x13ec;
c_in.power = 0x0;
c_in.control1 = 0x0;
#define RATE2000
#ifdef RATE2000
c_in.rate[0] = 0x0028; /* 2000 bit/s, no FEC */
c_in.rate[1] = 0x0000;
c_in.rate[2] = 0x0000;
c_in.rate[3] = 0x0000;
c_in.rate[4] = 0x6248;
#endif
#ifdef RATE3600_1200
c_in.rate[0] = 0x5048; /* 3600 bit/s, 1200 bit/s FEC */
c_in.rate[1] = 0x0001;
c_in.rate[2] = 0x0000;
c_in.rate[3] = 0x2412;
c_in.rate[4] = 0x6860;
#endif
c_in.unused[0] = 0x0;
c_in.unused[1] = 0x0;
c_in.unused[2] = 0x0;
c_in.dtmf = 0x00ff;
c_in.control2 = 0x8000;
/* put codec in run mode */
write_dongle(fd, run_state_run, sizeof(run_state_run));
//write_dongle(fd, data_item_1, sizeof(data_item_1));
//write_dongle(fd, (char*)&c_in, sizeof(c_in));
state = MSGSTATE_HDR1;
header = msg_type = msg_length = n = length = 0;
c_msg = uc_msg = 0;
for(i=0; i<100000; i++) {
/*
We can only reliably read one byte at a time. Until I
realised this there was "much wailing and gnashing of
teeth". Trying to read() n bytes read() returns n but may
actually reads some number between 1 and n. So it may only
read 1 byte int data[] but return n.
*/
r = read(fd, &data, 1);
assert(r == 1);
/* used state machine design from ambetest103.zip, SerialPort.cpp */
next_state = state;
switch(state) {
case MSGSTATE_HDR1:
header = data;
next_state = MSGSTATE_HDR2;
break;
case MSGSTATE_HDR2:
header |= data<<8;
msg_length = header & LENGTH_MASK;
msg_type = header & TYPE_MASK;
//printf("%0x %d\n", msg_type, msg_length);
if (length == 2) {
parse_message(msg_type, msg_length, msg_data);
next_state = MSGSTATE_HDR1;
}
else {
if (msg_length == 0x0)
length = 8192;
else
length = msg_length - 2;
n = 0;
next_state = MSGSTATE_DATA;
}
break;
case MSGSTATE_DATA:
msg_data[n++] = data;
length--;
if (length == 0) {
parse_message(msg_type, msg_length, msg_data);
next_state = MSGSTATE_HDR1;
}
break;
}
state = next_state;
}
printf("finished, c_msg = %d uc_msg = %d\n", c_msg, uc_msg);
write_dongle(fd, run_state_stop, sizeof(run_state_stop));
close(fd);
if (fin != NULL)
fclose(fin);
if (fout != NULL)
fclose(fout);
fclose(f);
return 0;
}

299
libs/libcodec2/unittest/ge_train.c Normal file
View File

@ -0,0 +1,299 @@
/*
ge_train.c
Jean Marc Valin Feb 2012
Joint pitch and energy VQ training program
usage:
cat GE | ./ge_train 2 1000000 8 > quantized
The first column is the log2 of the pitch compared to the lowest freq,
so log2(wo/pi*4000/50) where wo is the frequency your patch outputs. The
second column is the energy in dB, so 10*log10(1e-4+E)
*/
/*
Copyright (C) 2012 Jean-Marc Valin
All rights reserved.
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License version 2, as
published by the Free Software Foundation. This program is
distributed in the hope that it will be useful, but WITHOUT ANY
WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
License for more details.
You should have received a copy of the GNU Lesser General Public License
along with this program; if not, see <http://www.gnu.org/licenses/>.
*/
#include <valgrind/memcheck.h>
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#define MIN(a,b) ((a)<(b)?(a):(b))
//#define COEF 0.0
static float COEF[2] = {0.8, 0.9};
//static float COEF[2] = {0.0, 0.};
#define MAX_ENTRIES 16384
void compute_weights2(const float *x, const float *xp, float *w, int ndim)
{
w[0] = 30;
w[1] = 1;
if (x[1]<0)
{
w[0] *= .6;
w[1] *= .3;
}
if (x[1]<-10)
{
w[0] *= .3;
w[1] *= .3;
}
/* Higher weight if pitch is stable */
if (fabs(x[0]-xp[0])<.2)
{
w[0] *= 2;
w[1] *= 1.5;
} else if (fabs(x[0]-xp[0])>.5) /* Lower if not stable */
{
w[0] *= .5;
}
/* Lower weight for low energy */
if (x[1] < xp[1]-10)
{
w[1] *= .5;
}
if (x[1] < xp[1]-20)
{
w[1] *= .5;
}
//w[0] = 30;
//w[1] = 1;
/* Square the weights because it's applied on the squared error */
w[0] *= w[0];
w[1] *= w[1];
}
int find_nearest_weighted(const float *codebook, int nb_entries, float *x, const float *w, int ndim)
{
int i, j;
float min_dist = 1e15;
int nearest = 0;
for (i=0;i<nb_entries;i++)
{
float dist=0;
for (j=0;j<ndim;j++)
dist += w[j]*(x[j]-codebook[i*ndim+j])*(x[j]-codebook[i*ndim+j]);
if (dist<min_dist)
{
min_dist = dist;
nearest = i;
}
}
return nearest;
}
int quantize_ge(const float *x, const float *codebook1, int nb_entries, float *xq, int ndim)
{
int i, n1;
float err[ndim];
float w[ndim];
compute_weights2(x, xq, w, ndim);
for (i=0;i<ndim;i++)
err[i] = x[i]-COEF[i]*xq[i];
n1 = find_nearest_weighted(codebook1, nb_entries, err, w, ndim);
for (i=0;i<ndim;i++)
{
xq[i] = COEF[i]*xq[i] + codebook1[ndim*n1+i];
err[i] -= codebook1[ndim*n1+i];
}
return 0;
}
void split(float *codebook, int nb_entries, int ndim)
{
int i,j;
for (i=0;i<nb_entries;i++)
{
for (j=0;j<ndim;j++)
{
float delta = .01*(rand()/(float)RAND_MAX-.5);
codebook[i*ndim+j] += delta;
codebook[(i+nb_entries)*ndim+j] = codebook[i*ndim+j] - delta;
}
}
}
void update_weighted(float *data, float *weight, int nb_vectors, float *codebook, int nb_entries, int ndim)
{
int i,j;
float count[MAX_ENTRIES][ndim];
int nearest[nb_vectors];
//fprintf(stderr, "weighted: %d %d\n", nb_entries, ndim);
for (i=0;i<nb_entries;i++)
for (j=0;j<ndim;j++)
count[i][j] = 0;
for (i=0;i<nb_vectors;i++)
{
nearest[i] = find_nearest_weighted(codebook, nb_entries, data+i*ndim, weight+i*ndim, ndim);
}
for (i=0;i<nb_entries*ndim;i++)
codebook[i] = 0;
for (i=0;i<nb_vectors;i++)
{
int n = nearest[i];
for (j=0;j<ndim;j++)
{
float w = sqrt(weight[i*ndim+j]);
count[n][j]+=w;
codebook[n*ndim+j] += w*data[i*ndim+j];
}
}
//float w2=0;
for (i=0;i<nb_entries;i++)
{
for (j=0;j<ndim;j++)
codebook[i*ndim+j] *= (1./count[i][j]);
//w2 += (count[i]/(float)nb_vectors)*(count[i]/(float)nb_vectors);
}
//fprintf(stderr, "%f / %d\n", 1./w2, nb_entries);
}
void vq_train_weighted(float *data, float *weight, int nb_vectors, float *codebook, int nb_entries, int ndim)
{
int i, j, e;
e = 1;
for (j=0;j<ndim;j++)
codebook[j] = 0;
for (i=0;i<nb_vectors;i++)
for (j=0;j<ndim;j++)
codebook[j] += data[i*ndim+j];
for (j=0;j<ndim;j++)
codebook[j] *= (1./nb_vectors);
while (e< nb_entries)
{
#if 1
split(codebook, e, ndim);
e<<=1;
#else
split1(codebook, e, data, nb_vectors, ndim);
e++;
#endif
fprintf(stderr, "%d\n", e);
for (j=0;j<10;j++)
update_weighted(data, weight, nb_vectors, codebook, e, ndim);
}
}
int main(int argc, char **argv)
{
int i,j;
int nb_vectors, nb_entries, ndim;
float *data, *pred, *codebook, *codebook2, *codebook3;
float *weight, *weight2, *weight3;
float *delta;
double err[2] = {0, 0};
double werr[2] = {0, 0};
double wsum[2] = {0, 0};
ndim = atoi(argv[1]);
nb_vectors = atoi(argv[2]);
nb_entries = 1<<atoi(argv[3]);
data = malloc(nb_vectors*ndim*sizeof(*data));
weight = malloc(nb_vectors*ndim*sizeof(*weight));
weight2 = malloc(nb_vectors*ndim*sizeof(*weight2));
weight3 = malloc(nb_vectors*ndim*sizeof(*weight3));
pred = malloc(nb_vectors*ndim*sizeof(*pred));
codebook = malloc(nb_entries*ndim*sizeof(*codebook));
codebook2 = malloc(nb_entries*ndim*sizeof(*codebook2));
codebook3 = malloc(nb_entries*ndim*sizeof(*codebook3));
for (i=0;i<nb_vectors;i++)
{
if (feof(stdin))
break;
for (j=0;j<ndim;j++)
{
scanf("%f ", &data[i*ndim+j]);
}
}
nb_vectors = i;
VALGRIND_CHECK_MEM_IS_DEFINED(data, nb_entries*ndim);
for (i=0;i<nb_vectors;i++)
{
if (i==0)
compute_weights2(data+i*ndim, data+i*ndim, weight+i*ndim, ndim);
else
compute_weights2(data+i*ndim, data+(i-1)*ndim, weight+i*ndim, ndim);
}
for (i=0;i<ndim;i++)
pred[i] = data[i];
for (i=1;i<nb_vectors;i++)
{
for (j=0;j<ndim;j++)
pred[i*ndim+j] = data[i*ndim+j] - COEF[j]*data[(i-1)*ndim+j];
}
VALGRIND_CHECK_MEM_IS_DEFINED(pred, nb_entries*ndim);
vq_train_weighted(pred, weight, nb_vectors, codebook, nb_entries, ndim);
printf("%d %d\n", ndim, nb_entries);
for (i=0;i<nb_entries;i++)
{
for (j=0;j<ndim;j++)
{
printf("%f ", codebook[i*ndim+j]);
}
printf("\n");
}
delta = malloc(nb_vectors*ndim*sizeof(*data));
float xq[2] = {0,0};
for (i=0;i<nb_vectors;i++)
{
//int nearest = find_nearest_weighted(codebook, nb_entries, &pred[i*ndim], &weight[i*ndim], ndim);
quantize_ge(&data[i*ndim], codebook, nb_entries, xq, ndim);
//printf("%f %f\n", xq[0], xq[1]);
for (j=0;j<ndim;j++)
{
delta[i*ndim+j] = xq[j]-data[i*ndim+j];
err[j] += (delta[i*ndim+j])*(delta[i*ndim+j]);
werr[j] += weight[i*ndim+j]*(delta[i*ndim+j])*(delta[i*ndim+j]);
wsum[j] += weight[i*ndim+j];
//delta[i*ndim+j] = pred[i*ndim+j] - codebook[nearest*ndim+j];
//printf("%f ", delta[i*ndim+j]);
//err[j] += (delta[i*ndim+j])*(delta[i*ndim+j]);
}
//printf("\n");
}
fprintf(stderr, "GE RMS error: %f %f\n", sqrt(err[0]/nb_vectors), sqrt(err[1]/nb_vectors));
fprintf(stderr, "Weighted GE error: %f %f\n", sqrt(werr[0]/wsum[0]), sqrt(werr[1]/wsum[1]));
return 0;
}

BIN
libs/libcodec2/unittest/genampdata Executable file
View File

Binary file not shown.

104
libs/libcodec2/unittest/genampdata.c Normal file
View File

@ -0,0 +1,104 @@
/*
genampdata.c
Generates test sparse amplitude data for vqtrainsp testing.
*/
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <ctype.h>
#include <assert.h>
#include "../src/defines.h"
typedef struct {
float real;
float imag;
} COMP;
#define NVEC 200000
#define D 2
#define E 8
int main(void) {
FILE *f=fopen("testamp.flt", "wb");
int i, j, m, L, index;
float amp, noisey_amp, pitch, Wo;
float sparse_pe[MAX_AMP];
#ifdef TEST1
/* D fixed amplitude vectors of E elements long,
with D=2, E=8:
$ ./vqtrainsp testamp.flt 2 8 test.txt
test.txt should be same as training data.
*/
for(i=0; i<E; i++) {
amp = i+1;
for(j=0; j<D; j++)
fwrite(&amp, sizeof(float), 1, f);
}
#endif
#ifdef TEST2
/*
Bunch of amps uniformly distributed between -1 and 1. With e
entry "codebook" (1 dimensional vector or scalar):
$ ./vqtrainsp testamp.flt 1 e test.txt
should get std dev of 1/(e*sqrt(3))
*/
for(i=0; i<NVEC; i++) {
amp = 1.0 - 2.0*rand()/RAND_MAX;
fwrite(&amp, sizeof(float), 1, f);
}
#endif
#define TEST3
#ifdef TEST3
/*
Data for testing training of spare amplitudes. Similar to TEST1, each
sparse vector is set to the same amplitude.
/vqtrainsp testamp.flt 20 8 test.txt
*/
for(i=0; i<NVEC; i++) {
for(amp=1.0; amp<=8.0; amp++) {
pitch = P_MIN + (P_MAX-P_MIN)*((float)rand()/RAND_MAX);
Wo = TWO_PI/pitch;
L = floor(PI/Wo);
//printf("pitch %f Wo %f L %d\n", pitch, Wo, L);
for(m=0; m<MAX_AMP; m++) {
sparse_pe[m] = 0.0;
}
for(m=1; m<=L; m++) {
index = MAX_AMP*m*Wo/PI;
assert(index < MAX_AMP);
noisey_amp = amp + 0.2*(1.0 - 2.0*rand()/RAND_MAX);
sparse_pe[index] = noisey_amp;
#ifdef DBG
if (m < MAX_AMP/8)
printf(" %4.3f ", noisey_amp);
#endif
}
#ifdef DBG
printf("\n");
#endif
fwrite(sparse_pe, sizeof(float), MAX_AMP/8, f);
}
}
#endif
return 0;
}

BIN
libs/libcodec2/unittest/genphdata Executable file
View File

Binary file not shown.

93
libs/libcodec2/unittest/genphdata.c Normal file
View File

@ -0,0 +1,93 @@
/*
genphdata.c
Generates test phase data for vqtrainph testing.
*/
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <ctype.h>
#include <assert.h>
#include "../src/defines.h"
typedef struct {
float real;
float imag;
} COMP;
#define NVEC 100000
#define D 2
#define E 8
int main(void) {
FILE *f=fopen("testph.flt", "wb");
int i, m, L, index;
float angle, noisey_angle, pitch, Wo;
COMP c;
COMP sparse_pe[MAX_AMP];
#ifdef TEST1
for(i=0; i<D*E; i++) {
c.real = cos(i*TWO_PI/(M*D));
c.imag = sin(i*TWO_PI/(M*D));
fwrite(&c, sizeof(COMP), 1, f);
}
#endif
#ifdef TEST2
/*
Bunch of random phases, should get std dev per element of
pi/(sqrt(3)*pow(2,b/D)), or 0.321 for (b=5, D=2):
./vqtrainph testph.flt 2 32 test.txt
*/
for(i=0; i<NVEC; i++) {
angle = PI*(1.0 - 2.0*rand()/RAND_MAX);
c.real = cos(angle);
c.imag = sin(angle);
fwrite(&c, sizeof(COMP), 1, f);
}
#endif
#define TEST3
#ifdef TEST3
/*
Data for testing training in sparse phases. No correlation, so
should be same performance as TEST2. Attempting to train a
MAX_AMP/4 = 20 (first 1 kHz) phase quantiser.
*/
angle = 0;
for(i=0; i<NVEC; i++) {
pitch = P_MIN + (P_MAX-P_MIN)*((float)rand()/RAND_MAX);
//pitch = 40;
Wo = TWO_PI/pitch;
L = floor(PI/Wo);
//printf("pitch %f Wo %f L %d\n", pitch, Wo, L);
for(m=0; m<MAX_AMP; m++) {
sparse_pe[m].real = 0.0;
sparse_pe[m].imag = 0.0;
}
angle += PI/8;
for(m=1; m<=L; m++) {
noisey_angle = angle + (PI/16)*(1.0 - 2.0*rand()/RAND_MAX);
//angle = (PI/16)*(1.0 - 2.0*rand()/RAND_MAX);
index = MAX_AMP*m*Wo/PI;
assert(index < MAX_AMP);
sparse_pe[index].real = cos(noisey_angle);
sparse_pe[index].imag = sin(noisey_angle);
}
fwrite(&sparse_pe, sizeof(COMP), MAX_AMP/4, f);
}
#endif
return 0;
}

BIN
libs/libcodec2/unittest/lspsync Executable file
View File

Binary file not shown.

115
libs/libcodec2/unittest/lspsync.c Normal file
View File

@ -0,0 +1,115 @@
/*
lspsync.c
David Rowe 24 May 2012
Attempt at using LSP information to provide frame sync. If we have
correct frame alignment, LSPs will not need sorting.
However this method as tested appears unreliable, often several
sync positions per frame are found, even with a F=10 memory. For
F=6, about 87% relaible. This might be useful if combined with a
another sync method, for example a single alternating sync bit per
frame.
*/
#include <assert.h>
#include <stdio.h>
#include <stdlib.h>
#include "codec2.h"
#include "defines.h"
#include "quantise.h"
#define F 6 /* look at LSP ordering in F-1 frames */
#define CORRECT_OFFSET 10 /* LSPs start 10 bits int frame qt 2400 bit/s */
static int check_candidate(char bits[], int offset)
{
int i;
int lsp_indexes[LPC_ORD];
float lsps[LPC_ORD];
unsigned int nbit = offset;
int swaps;
for(i=0; i<LSP_SCALAR_INDEXES; i++) {
lsp_indexes[i] = unpack(bits, &nbit, lsp_bits(i));
}
decode_lsps_scalar(lsps, lsp_indexes, LPC_ORD);
swaps = check_lsp_order(lsps, LPC_ORD);
return swaps;
}
int main(int argc, char *argv[]) {
struct CODEC2 *c2;
int i,offset, nsamples, nbits, nbytes, frames;
short *speech;
char *bits;
FILE *fin;
int swaps, pass, fail, match;
c2 = codec2_create(CODEC2_MODE_2400);
nsamples = codec2_samples_per_frame(c2);
nbits = codec2_bits_per_frame(c2);
nbytes = nbits/8;
speech = (short*)malloc(nsamples*sizeof(short));
/* keep FRAMES frame memory of bit stream */
bits = (char*)malloc(F*nbytes*sizeof(char));
for(i=0; i<F*nbytes; i++)
bits[i] = 0;
fin = fopen("../raw/hts1a.raw", "rb");
assert(fin != NULL);
match = pass = fail = frames = 0;
/* prime memeory with first frame to ensure we don't start
checking until we have two frames of coded bits */
fread(speech, sizeof(short), nsamples, fin);
frames++;
codec2_encode(c2, &bits[(F-2)*nbytes], speech);
/* OK start looking for correct frame offset */
while(fread(speech, sizeof(short), nsamples, fin) == nsamples) {
frames++;
codec2_encode(c2, &bits[(F-1)*nbytes], speech);
for(offset=0; offset<nbits; offset++) {
swaps = check_candidate(bits, offset);
if (swaps == 0) {
/* OK found a candidate .. lets check a F-1 frames in total */
for(i=0; i<(F-1); i++)
swaps += check_candidate(bits, offset + nbits*i);
if (swaps == 0) {
printf("frame %d offset: %d swaps: %d\n", frames, offset, swaps);
match++;
if (offset == CORRECT_OFFSET)
pass++;
else
fail++;
}
}
}
/* update F frame memory of bits */
for(i=0; i<nbytes*(F-1); i++)
bits[i] = bits[i+nbytes];
}
fclose(fin);
free(speech);
free(bits);
codec2_destroy(c2);
printf("passed %f %%\n", (float)pass*100.0/match);
return 0;
}

38
libs/libcodec2/unittest/mksine.c Normal file
View File

@ -0,0 +1,38 @@
/*
mksine.c
David Rowe
10 Nov 2010
Creates a file of sine wave samples.
*/
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#define TWO_PI 6.283185307
#define N 8000
#define FS 8000.0
#define AMP 1000.0
int main(int argc, char *argv[]) {
FILE *f;
int i;
float freq;
short buf[N];
if (argc != 3) {
printf("usage: %s outputFile frequencyHz\n", argv[0]);
exit(1);
}
f = fopen(argv[1] ,"wb");
freq = atof(argv[2]);
for(i=0; i<N; i++)
buf[i] = AMP*cos(freq*i*(TWO_PI/FS));
fwrite(buf, sizeof(short), N, f);
return 0;
}

BIN
libs/libcodec2/unittest/polar2rect Executable file
View File

Binary file not shown.

54
libs/libcodec2/unittest/polar2rect.c Normal file
View File

@ -0,0 +1,54 @@
/*
polar2rect.c
David Rowe 28 July 2013
Convert a file of sparse phases in polar (angle) format to a file in rect
format.
*/
#include <assert.h>
#include <math.h>
#include <stdlib.h>
#include <stdio.h>
typedef struct {
float real;
float imag;
} COMP;
int main(int argc, char *argv[]) {
FILE *fpolar;
FILE *frect;
float polar;
COMP rect;
if (argc != 3) {
printf("usage: %s polarFile rectFile\n", argv[0]);
exit(0);
}
fpolar = fopen(argv[1], "rb");
assert(fpolar != NULL);
frect = fopen(argv[2], "wb");
assert(frect != NULL);
while (fread(&polar, sizeof(float), 1, fpolar) != 0) {
if (polar == 0.0) {
/* this values indicates the VQ training should ignore
this vector element. It's not a valid phase as it
doesn't have mangitude of 1.0 */
rect.real = 0.0;
rect.imag = 0.0;
}
else {
rect.real = cos(polar);
rect.imag = sin(polar);
}
fwrite(&rect, sizeof(COMP), 1, frect);
}
fclose(fpolar);
fclose(frect);
return 0;
}

BIN
libs/libcodec2/unittest/pre Executable file
View File

Binary file not shown.

59
libs/libcodec2/unittest/pre.c Normal file
View File

@ -0,0 +1,59 @@
/*
pre.c
David Rowe
Sep 26 2012
Takes audio from a file, pre-emphasises, and sends to output file.
*/
#include <assert.h>
#include <math.h>
#include <errno.h>
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include "lpc.h"
#define N 80
int main(int argc, char*argv[]) {
FILE *fin, *fout;
short buf[N];
float Sn[N], Sn_pre[N];
float pre_mem = 0.0;
int i;
if (argc != 3) {
printf("usage: pre InputRawSpeechFile OutputRawSpeechFile\n");
printf("e.g pre input.raw output.raw");
exit(1);
}
if (strcmp(argv[1], "-") == 0) fin = stdin;
else if ( (fin = fopen(argv[1],"rb")) == NULL ) {
fprintf(stderr, "Error opening input speech file: %s: %s.\n",
argv[1], strerror(errno));
exit(1);
}
if (strcmp(argv[2], "-") == 0) fout = stdout;
else if ( (fout = fopen(argv[2],"wb")) == NULL ) {
fprintf(stderr, "Error opening output speech file: %s: %s.\n",
argv[2], strerror(errno));
exit(1);
}
while(fread(buf, sizeof(short), N, fin) == N) {
for(i=0; i<N; i++)
Sn[i] = buf[i];
pre_emp(Sn_pre, Sn, &pre_mem, N);
for(i=0; i<N; i++)
buf[i] = Sn_pre[i];
fwrite(buf, sizeof(short), N, fout);
}
fclose(fin);
fclose(fout);
return 0;
}

27
libs/libcodec2/unittest/run_tests.sh Normal file
View File

@ -0,0 +1,27 @@
#!/bin/bash
# verbose, skittish
set -ex
# run valgrind?
VG="valgrind --error-exitcode=149 --leak-check=full --show-reachable=yes"
# make sure we're up-to-date
make
# runs, not certaion wht it does yet
${VG} ./tcodec2
${VG} ./tinterp
${VG} ./tquant
# these fail, missing arguments
${VG} ./extract
${VG} ./genlsp
${VG} ./genres
${VG} ./scalarlsptest
${VG} ./tnlp
${VG} ./vq_train_jvm
${VG} ./vqtrain
${VG} ./vqtrainjnd

BIN
libs/libcodec2/unittest/scalarlsptest Executable file
View File

Binary file not shown.

106
libs/libcodec2/unittest/scalarlsptest.c Normal file
View File

@ -0,0 +1,106 @@
/*---------------------------------------------------------------------------*\
FILE........: scalarlsptest.c
AUTHOR......: David Rowe
DATE CREATED: 8/2/12
Test Scalar LSP quantiser, output variance of quantisation error.
\*---------------------------------------------------------------------------*/
/*
Copyright (C) 2012 David Rowe
All rights reserved.
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License version 2, as
published by the Free Software Foundation. This program is
distributed in the hope that it will be useful, but WITHOUT ANY
WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
License for more details.
You should have received a copy of the GNU Lesser General Public License
along with this program; if not, see <http://www.gnu.org/licenses/>.
*/
#include <assert.h>
#include <stdio.h>
#include <stdlib.h>
#include "defines.h"
#include "quantise.h"
/*---------------------------------------------------------------------------*\
MAIN
\*---------------------------------------------------------------------------*/
int main(int argc, char *argv[])
{
FILE *ftrain; /* LSP training data text file */
float lsp[LPC_ORD]; /* LSP input vector in rads */
float lsp_hz[LPC_ORD]; /* LSP input vector in Hz */
int vectors; /* number LSP vectors processed */
int k,m; /* LSP vector order and codebook size */
int index;
float wt[1]; /* weighting (not used here for scalars) */
const float *cb; /* LSP quantiser codebook */
int i, ret;
float total_se;
if (argc < 2) {
printf("usage: %s InputFile\n", argv[0]);
exit(1);
}
if ((ftrain = fopen(argv[1],"rt")) == NULL) {
printf("Error opening input file: %s\n",argv[1]);
exit(0);
}
total_se = 0.0;
vectors = 0;
wt[0] = 1.0;
/* Main loop */
while(!feof(ftrain)) {
/* Read LSP input vector speech */
for (i=0; i<LPC_ORD; i++) {
ret = fscanf(ftrain, "%f ", &lsp[i]);
}
vectors++;
if ((vectors % 1000) == 0)
printf("\r%d vectors", vectors);
/* convert from radians to Hz so we can use human readable
frequencies */
for(i=0; i<LPC_ORD; i++)
lsp_hz[i] = (4000.0/PI)*lsp[i];
/* simple uniform scalar quantisers */
for(i=0; i<LPC_ORD; i++) {
k = lsp_cb[i].k;
m = lsp_cb[i].m;
cb = lsp_cb[i].cb;
index = quantise(cb, &lsp_hz[i], wt, k, m, &total_se);
//printf("k %d m %d lsp[%d] %f %f se %f\n", k,m,i,lsp_hz[i], cb[index],se);
}
//printf("total se %f\n", total_se);
//exit(0);
}
fclose(ftrain);
printf("\n variance = %f\n", ((PI*PI)/(4000.0*4000.0))*total_se/vectors);
return 0;
}

176
libs/libcodec2/unittest/speexlsptest.c Normal file
View File

@ -0,0 +1,176 @@
/*---------------------------------------------------------------------------*\
FILE........: speexlsptest.c
AUTHOR......: David Rowe
DATE CREATED: 24/8/09
Test LPC to LSP conversion and quantisation using Speex LSP quantiser.
\*---------------------------------------------------------------------------*/
/*
Copyright (C) 2009 David Rowe
All rights reserved.
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License version 2, as
published by the Free Software Foundation. This program is
distributed in the hope that it will be useful, but WITHOUT ANY
WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
License for more details.
You should have received a copy of the GNU Lesser General Public License
along with this program; if not, see <http://www.gnu.org/licenses/>.
*/
#include <assert.h>
#include <stdio.h>
#include <stdlib.h>
#include <lpc.h>
#include <lsp.h>
#include <sd.h>
#define N 160
#define P 10
#define LPC_FLOOR 0.0002 /* autocorrelation floor */
#define LSP_DELTA1 0.2 /* grid spacing for LSP root searches */
#define NDFT 256 /* DFT size for SD calculation */
/* Speex lag window */
const float lag_window[11] = {
1.00000, 0.99716, 0.98869, 0.97474, 0.95554, 0.93140, 0.90273, 0.86998,
0.83367, 0.79434, 0.75258
};
/*---------------------------------------------------------------------------*\
find_aks_for_lsp()
This function takes a frame of samples, and determines the linear
prediction coefficients for that frame of samples. Modified version of
find_aks from lpc.c to include autocorrelation noise floor and lag window
to match Speex processing steps prior to LSP conversion.
\*---------------------------------------------------------------------------*/
void find_aks_for_lsp(
float Sn[], /* Nsam samples with order sample memory */
float a[], /* order+1 LPCs with first coeff 1.0 */
int Nsam, /* number of input speech samples */
int order, /* order of the LPC analysis */
float *E /* residual energy */
)
{
float Wn[N]; /* windowed frame of Nsam speech samples */
float R[P+1]; /* order+1 autocorrelation values of Sn[] */
int i;
hanning_window(Sn,Wn,Nsam);
autocorrelate(Wn,R,Nsam,order);
R[0] += LPC_FLOOR;
assert(order == 10); /* lag window only defined for order == 10 */
for(i=0; i<=order; i++)
R[i] *= lag_window[i];
levinson_durbin(R,a,order);
*E = 0.0;
for(i=0; i<=order; i++)
*E += a[i]*R[i];
if (*E < 0.0)
*E = 1E-12;
}
/*---------------------------------------------------------------------------*\
MAIN
\*---------------------------------------------------------------------------*/
int main(int argc, char *argv[])
{
FILE *fin; /* input speech files */
short buf[N]; /* buffer of 16 bit speech samples */
float Sn[P+N]; /* input speech samples */
float E;
float ak[P+1]; /* LP coeffs */
float ak_[P+1]; /* quantised LP coeffs */
float lsp[P];
float lsp_[P]; /* quantised LSPs */
int roots; /* number of LSP roots found */
int frames; /* frames processed so far */
int i; /* loop variables */
SpeexBits bits;
float sd; /* SD for this frame */
float totsd; /* accumulated SD so far */
int gt2,gt4; /* number of frames > 2 and 4 dB SD */
int unstables; /* number of unstable LSP frames */
if (argc < 2) {
printf("usage: %s InputFile\n", argv[0]);
exit(0);
}
/* Open files */
if ((fin = fopen(argv[1],"rb")) == NULL) {
printf("Error opening input file: %s\n",argv[1]);
exit(0);
}
/* Initialise */
frames = 0;
for(i=0; i<P; i++) {
Sn[i] = 0.0;
}
ak_[0] = 1.0;
speex_bits_init(&bits);
totsd = 0.0;
unstables = 0;
gt2 = 0; gt4 = 0;
/* Main loop */
while( (fread(buf,sizeof(short),N,fin)) == N) {
frames++;
for(i=0; i<N; i++)
Sn[P+i] = (float)buf[i];
/* convert to LSP domain and back */
find_aks(&Sn[P], ak, N, P, &E);
roots = lpc_to_lsp(&ak[1], P , lsp, 10, LSP_DELTA1, NULL);
if (roots == P) {
speex_bits_reset(&bits);
lsp_quant_lbr(lsp, lsp_, P, &bits);
lsp_to_lpc(lsp_, &ak_[1], P, NULL);
/* measure spectral distortion */
sd = spectral_dist(ak, ak_, P, NDFT);
if (sd > 2.0) gt2++;
if (sd > 4.0) gt4++;
totsd += sd;
}
else
unstables++;
}
fclose(fin);
printf("frames = %d Av sd = %3.2f dB", frames, totsd/frames);
printf(" >2 dB %3.2f%% >4 dB %3.2f%% unstables: %d\n",gt2*100.0/frames,
gt4*100.0/frames, unstables);
return 0;
}

BIN
libs/libcodec2/unittest/t48_8 Executable file
View File

Binary file not shown.

112
libs/libcodec2/unittest/t48_8.c Normal file
View File

@ -0,0 +1,112 @@
/*
t48_8.c
David Rowe
May 10 2012
Unit test for 48 to 8 kHz sample rate conversion functions. I
evaluated output by plotting using Octave and looking for jaggies:
pl("../unittest/out48.raw",1,3000)
pl("../unittest/out8.raw",1,3000)
Listening to it also shows up anything nasty:
$ play -s -2 -r 48000 out48.raw
$ play -s -2 -r 8000 out8.raw
*/
#include <assert.h>
#include <math.h>
#include <stdlib.h>
#include <stdio.h>
#include "fdmdv.h"
#define N8 160 /* procssing buffer size at 8 kHz */
#define N48 (N8*FDMDV_OS)
#define MEM8 (FDMDV_OS_TAPS/FDMDV_OS)
#define FRAMES 50
#define TWO_PI 6.283185307
#define FS 8000
#define SINE
int main() {
float in8k[MEM8 + N8];
float out48k[N48];
short out48k_short[N48];
FILE *f48;
float in48k[FDMDV_OS_TAPS + N48];
float out8k[N48];
short out8k_short[N8];
FILE *f8;
int i,f,t,t1;
float freq = 800.0;
f48 = fopen("out48.raw", "wb");
assert(f48 != NULL);
f8 = fopen("out8.raw", "wb");
assert(f8 != NULL);
/* clear filter memories */
for(i=0; i<MEM8; i++)
in8k[i] = 0.0;
for(i=0; i<FDMDV_OS_TAPS; i++)
in48k[i] = 0.0;
t = t1 = 0;
for(f=0; f<FRAMES; f++) {
#ifdef DC
for(i=0; i<N8; i++)
in8k[MEM8+i] = 16000.0;
#endif
#ifdef SINE
for(i=0; i<N8; i++,t++)
in8k[MEM8+i] = 16000.0*cos(TWO_PI*t*freq/FS);
#endif
/* upsample */
fdmdv_8_to_48(out48k, &in8k[MEM8], N8);
/*
for(i=0; i<MEM8; i++)
in8k[i] = in8k[i+N8];
*/
/* save 48k to disk for plotting and check out */
for(i=0; i<N48; i++)
out48k_short[i] = (short)out48k[i];
fwrite(out48k_short, sizeof(short), N48, f48);
/* add a 10 kHz spurious signal for fun, we want down sampler to
knock this out */
for(i=0; i<N48; i++,t1++)
in48k[i+FDMDV_OS_TAPS] = out48k[i] + 16000.0*cos(TWO_PI*t1*1E4/FS);
/* downsample */
fdmdv_48_to_8(out8k, &in48k[FDMDV_OS_TAPS], N8);
/*
for(i=0; i<FDMDV_OS_TAPS; i++)
in48k[i] = in48k[i+N48];
*/
/* save 8k to disk for plotting and check out */
for(i=0; i<N8; i++)
out8k_short[i] = (short)out8k[i];
fwrite(out8k_short, sizeof(short), N8, f8);
}
fclose(f48);
fclose(f8);
return 0;
}

BIN
libs/libcodec2/unittest/tfdmdv Executable file
View File

Binary file not shown.

282
libs/libcodec2/unittest/tfdmdv.c Normal file
View File

@ -0,0 +1,282 @@
/*---------------------------------------------------------------------------*\
FILE........: tfdmdv.c
AUTHOR......: David Rowe
DATE CREATED: April 16 2012
Tests for the C version of the FDMDV modem. This program outputs a
file of Octave vectors that are loaded and automatically tested
against the Octave version of the modem by the Octave script
tfmddv.m
\*---------------------------------------------------------------------------*/
/*
Copyright (C) 2012 David Rowe
All rights reserved.
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License version 2, as
published by the Free Software Foundation. This program is
distributed in the hope that it will be useful, but WITHOUT ANY
WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
License for more details.
You should have received a copy of the GNU Lesser General Public License
along with this program; if not, see <http://www.gnu.org/licenses/>.
*/
#include <assert.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include "fdmdv_internal.h"
#include "fdmdv.h"
#include "octave.h"
#define FRAMES 25
#define CHANNEL_BUF_SIZE (10*M)
int main(int argc, char *argv[])
{
struct FDMDV *fdmdv;
int tx_bits[FDMDV_BITS_PER_FRAME];
COMP tx_symbols[NC+1];
COMP tx_baseband[NC+1][M];
COMP tx_fdm[M];
float channel[CHANNEL_BUF_SIZE];
int channel_count;
COMP rx_fdm[M+M/P];
float foff_coarse;
int nin, next_nin;
COMP rx_fdm_fcorr[M+M/P];
COMP rx_baseband[NC+1][M+M/P];
COMP rx_filt[NC+1][P+1];
float rx_timing;
float env[NT*P];
COMP rx_symbols[NC+1];
int rx_bits[FDMDV_BITS_PER_FRAME];
float foff_fine;
int sync_bit;
int fest_state;
int tx_bits_log[FDMDV_BITS_PER_FRAME*FRAMES];
COMP tx_symbols_log[(NC+1)*FRAMES];
COMP tx_baseband_log[(NC+1)][M*FRAMES];
COMP tx_fdm_log[M*FRAMES];
COMP pilot_baseband1_log[NPILOTBASEBAND*FRAMES];
COMP pilot_baseband2_log[NPILOTBASEBAND*FRAMES];
COMP pilot_lpf1_log[NPILOTLPF*FRAMES];
COMP pilot_lpf2_log[NPILOTLPF*FRAMES];
COMP S1_log[MPILOTFFT*FRAMES];
COMP S2_log[MPILOTFFT*FRAMES];
float foff_coarse_log[FRAMES];
float foff_log[FRAMES];
COMP rx_baseband_log[(NC+1)][(M+M/P)*FRAMES];
int rx_baseband_log_col_index;
COMP rx_filt_log[NC+1][(P+1)*FRAMES];
int rx_filt_log_col_index;
float env_log[NT*P*FRAMES];
float rx_timing_log[FRAMES];
COMP rx_symbols_log[NC+1][FRAMES];
float sig_est_log[NC+1][FRAMES];
float noise_est_log[NC+1][FRAMES];
int rx_bits_log[FDMDV_BITS_PER_FRAME*FRAMES];
float foff_fine_log[FRAMES];
int sync_bit_log[FRAMES];
int coarse_fine_log[FRAMES];
int nin_log[FRAMES];
FILE *fout;
int f,c,i,j;
fdmdv = fdmdv_create();
next_nin = M;
channel_count = 0;
rx_baseband_log_col_index = 0;
rx_filt_log_col_index = 0;
printf("sizeof FDMDV states: %d bytes\n", sizeof(struct FDMDV));
for(f=0; f<FRAMES; f++) {
/* --------------------------------------------------------*\
Modulator
\*---------------------------------------------------------*/
fdmdv_get_test_bits(fdmdv, tx_bits);
bits_to_dqpsk_symbols(tx_symbols, fdmdv->prev_tx_symbols, tx_bits, &fdmdv->tx_pilot_bit);
memcpy(fdmdv->prev_tx_symbols, tx_symbols, sizeof(COMP)*(NC+1));
tx_filter(tx_baseband, tx_symbols, fdmdv->tx_filter_memory);
fdm_upconvert(tx_fdm, tx_baseband, fdmdv->phase_tx, fdmdv->freq);
/* --------------------------------------------------------*\
Channel
\*---------------------------------------------------------*/
nin = next_nin;
/*
if (f == 2)
nin = 120;
if (f == 3)
nin = 200;
if ((f !=2) && (f != 3))
nin = M;
*/
/* add M tx samples to end of buffer */
assert((channel_count + M) < CHANNEL_BUF_SIZE);
for(i=0; i<M; i++)
channel[channel_count+i] = tx_fdm[i].real;
channel_count += M;
/* take nin samples from start of buffer */
for(i=0; i<nin; i++) {
rx_fdm[i].real = channel[i];
rx_fdm[i].imag = 0;
}
/* shift buffer back */
for(i=0,j=nin; j<channel_count; i++,j++)
channel[i] = channel[j];
channel_count -= nin;
/* --------------------------------------------------------*\
Demodulator
\*---------------------------------------------------------*/
/* freq offset estimation and correction */
foff_coarse = rx_est_freq_offset(fdmdv, rx_fdm, nin);
if (fdmdv->coarse_fine == COARSE)
fdmdv->foff = foff_coarse;
fdmdv_freq_shift(rx_fdm_fcorr, rx_fdm, fdmdv->foff, &fdmdv->foff_rect, &fdmdv->foff_phase_rect, nin);
/* baseband processing */
fdm_downconvert(rx_baseband, rx_fdm_fcorr, fdmdv->phase_rx, fdmdv->freq, nin);
rx_filter(rx_filt, rx_baseband, fdmdv->rx_filter_memory, nin);
rx_timing = rx_est_timing(rx_symbols, rx_filt, rx_baseband, fdmdv->rx_filter_mem_timing, env, fdmdv->rx_baseband_mem_timing, nin);
foff_fine = qpsk_to_bits(rx_bits, &sync_bit, fdmdv->phase_difference, fdmdv->prev_rx_symbols, rx_symbols);
snr_update(fdmdv->sig_est, fdmdv->noise_est, fdmdv->phase_difference);
memcpy(fdmdv->prev_rx_symbols, rx_symbols, sizeof(COMP)*(NC+1));
next_nin = M;
if (rx_timing > 2*M/P)
next_nin += M/P;
if (rx_timing < 0)
next_nin -= M/P;
fdmdv->coarse_fine = freq_state(sync_bit, &fdmdv->fest_state);
fdmdv->foff -= TRACK_COEFF*foff_fine;
/* --------------------------------------------------------*\
Log each vector
\*---------------------------------------------------------*/
memcpy(&tx_bits_log[FDMDV_BITS_PER_FRAME*f], tx_bits, sizeof(int)*FDMDV_BITS_PER_FRAME);
memcpy(&tx_symbols_log[(NC+1)*f], tx_symbols, sizeof(COMP)*(NC+1));
for(c=0; c<NC+1; c++)
for(i=0; i<M; i++)
tx_baseband_log[c][f*M+i] = tx_baseband[c][i];
memcpy(&tx_fdm_log[M*f], tx_fdm, sizeof(COMP)*M);
/* freq offset estimation */
memcpy(&pilot_baseband1_log[f*NPILOTBASEBAND], fdmdv->pilot_baseband1, sizeof(COMP)*NPILOTBASEBAND);
memcpy(&pilot_baseband2_log[f*NPILOTBASEBAND], fdmdv->pilot_baseband2, sizeof(COMP)*NPILOTBASEBAND);
memcpy(&pilot_lpf1_log[f*NPILOTLPF], fdmdv->pilot_lpf1, sizeof(COMP)*NPILOTLPF);
memcpy(&pilot_lpf2_log[f*NPILOTLPF], fdmdv->pilot_lpf2, sizeof(COMP)*NPILOTLPF);
memcpy(&S1_log[f*MPILOTFFT], fdmdv->S1, sizeof(COMP)*MPILOTFFT);
memcpy(&S2_log[f*MPILOTFFT], fdmdv->S2, sizeof(COMP)*MPILOTFFT);
foff_coarse_log[f] = foff_coarse;
foff_log[f] = fdmdv->foff;
/* rx down conversion */
for(c=0; c<NC+1; c++) {
for(i=0; i<nin; i++)
rx_baseband_log[c][rx_baseband_log_col_index + i] = rx_baseband[c][i];
}
rx_baseband_log_col_index += nin;
/* rx filtering */
for(c=0; c<NC+1; c++) {
for(i=0; i<(P*nin)/M; i++)
rx_filt_log[c][rx_filt_log_col_index + i] = rx_filt[c][i];
}
rx_filt_log_col_index += (P*nin)/M;
/* timing estimation */
memcpy(&env_log[NT*P*f], env, sizeof(float)*NT*P);
rx_timing_log[f] = rx_timing;
nin_log[f] = nin;
for(c=0; c<NC+1; c++)
rx_symbols_log[c][f] = rx_symbols[c];
/* qpsk_to_bits() */
memcpy(&rx_bits_log[FDMDV_BITS_PER_FRAME*f], rx_bits, sizeof(int)*FDMDV_BITS_PER_FRAME);
for(c=0; c<NC+1; c++) {
sig_est_log[c][f] = fdmdv->sig_est[c];
noise_est_log[c][f] = fdmdv->noise_est[c];
}
foff_fine_log[f] = foff_fine;
sync_bit_log[f] = sync_bit;
coarse_fine_log[f] = fdmdv->coarse_fine;
}
/*---------------------------------------------------------*\
Dump logs to Octave file for evaluation
by tfdmdv.m Octave script
\*---------------------------------------------------------*/
fout = fopen("tfdmdv_out.txt","wt");
assert(fout != NULL);
fprintf(fout, "# Created by tfdmdv.c\n");
octave_save_int(fout, "tx_bits_log_c", tx_bits_log, 1, FDMDV_BITS_PER_FRAME*FRAMES);
octave_save_complex(fout, "tx_symbols_log_c", tx_symbols_log, 1, (NC+1)*FRAMES, (NC+1)*FRAMES);
octave_save_complex(fout, "tx_baseband_log_c", (COMP*)tx_baseband_log, (NC+1), M*FRAMES, M*FRAMES);
octave_save_complex(fout, "tx_fdm_log_c", (COMP*)tx_fdm_log, 1, M*FRAMES, M*FRAMES);
octave_save_complex(fout, "pilot_lut_c", (COMP*)fdmdv->pilot_lut, 1, NPILOT_LUT, NPILOT_LUT);
octave_save_complex(fout, "pilot_baseband1_log_c", pilot_baseband1_log, 1, NPILOTBASEBAND*FRAMES, NPILOTBASEBAND*FRAMES);
octave_save_complex(fout, "pilot_baseband2_log_c", pilot_baseband2_log, 1, NPILOTBASEBAND*FRAMES, NPILOTBASEBAND*FRAMES);
octave_save_complex(fout, "pilot_lpf1_log_c", pilot_lpf1_log, 1, NPILOTLPF*FRAMES, NPILOTLPF*FRAMES);
octave_save_complex(fout, "pilot_lpf2_log_c", pilot_lpf2_log, 1, NPILOTLPF*FRAMES, NPILOTLPF*FRAMES);
octave_save_complex(fout, "S1_log_c", S1_log, 1, MPILOTFFT*FRAMES, MPILOTFFT*FRAMES);
octave_save_complex(fout, "S2_log_c", S2_log, 1, MPILOTFFT*FRAMES, MPILOTFFT*FRAMES);
octave_save_float(fout, "foff_log_c", foff_log, 1, FRAMES, FRAMES);
octave_save_float(fout, "foff_coarse_log_c", foff_coarse_log, 1, FRAMES, FRAMES);
octave_save_complex(fout, "rx_baseband_log_c", (COMP*)rx_baseband_log, (NC+1), rx_baseband_log_col_index, (M+M/P)*FRAMES);
octave_save_complex(fout, "rx_filt_log_c", (COMP*)rx_filt_log, (NC+1), rx_filt_log_col_index, (P+1)*FRAMES);
octave_save_float(fout, "env_log_c", env_log, 1, NT*P*FRAMES, NT*P*FRAMES);
octave_save_float(fout, "rx_timing_log_c", rx_timing_log, 1, FRAMES, FRAMES);
octave_save_complex(fout, "rx_symbols_log_c", (COMP*)rx_symbols_log, (NC+1), FRAMES, FRAMES);
octave_save_float(fout, "sig_est_log_c", (float*)sig_est_log, (NC+1), FRAMES, FRAMES);
octave_save_float(fout, "noise_est_log_c", (float*)noise_est_log, (NC+1), FRAMES, FRAMES);
octave_save_int(fout, "rx_bits_log_c", rx_bits_log, 1, FDMDV_BITS_PER_FRAME*FRAMES);
octave_save_float(fout, "foff_fine_log_c", foff_fine_log, 1, FRAMES, FRAMES);
octave_save_int(fout, "sync_bit_log_c", sync_bit_log, 1, FRAMES);
octave_save_int(fout, "coarse_fine_log_c", coarse_fine_log, 1, FRAMES);
octave_save_int(fout, "nin_log_c", nin_log, 1, FRAMES);
fclose(fout);
fdmdv_destroy(fdmdv);
return 0;
}

BIN
libs/libcodec2/unittest/tfifo Executable file
View File

Binary file not shown.

103
libs/libcodec2/unittest/tfifo.c Normal file
View File

@ -0,0 +1,103 @@
/*
tfifo.c
David Rowe
Nov 19 2012
Takes FIFOs, in particular thread safety.
*/
#include <assert.h>
#include <stdio.h>
#include <pthread.h>
#include "fifo.h"
#define FIFO_SZ 1024
#define WRITE_SZ 10
#define READ_SZ 8
#define N_MAX 100
#define LOOPS 1000000
int run_thread = 1;
struct FIFO *f;
void writer(void);
void *writer_thread(void *data);
pthread_mutex_t mutex = PTHREAD_MUTEX_INITIALIZER;
#define USE_THREADS
//#define USE_MUTEX
int main() {
pthread_t awriter_thread;
int i,j;
short read_buf[READ_SZ];
int n_out = 0;
int sucess;
f = fifo_create(FIFO_SZ);
#ifdef USE_THREADS
pthread_create(&awriter_thread, NULL, writer_thread, NULL);
#endif
for(i=0; i<LOOPS; ) {
#ifndef USE_THREADS
writer();
#endif
#ifdef USE_MUTEX
pthread_mutex_lock(&mutex);
#endif
sucess = (fifo_read(f, read_buf, READ_SZ) == 0);
#ifdef USE_MUTEX
pthread_mutex_unlock(&mutex);
#endif
if (sucess) {
for(j=0; j<READ_SZ; j++) {
if (read_buf[j] != n_out)
printf("error: %d %d\n", read_buf[j], n_out);
n_out++;
if (n_out == N_MAX)
n_out = 0;
}
i++;
}
}
#ifdef USE_THREADS
run_thread = 0;
pthread_join(awriter_thread,NULL);
#endif
return 0;
}
int n_in = 0;
void writer(void) {
short write_buf[WRITE_SZ];
int i;
if ((FIFO_SZ - fifo_used(f)) > WRITE_SZ) {
for(i=0; i<WRITE_SZ; i++) {
write_buf[i] = n_in++;
if (n_in == N_MAX)
n_in = 0;
}
#ifdef USE_MUTEX
pthread_mutex_lock(&mutex);
#endif
fifo_write(f, write_buf, WRITE_SZ);
pthread_mutex_unlock(&mutex);
}
}
void *writer_thread(void *data) {
while(run_thread) {
writer();
}
return NULL;
}

BIN
libs/libcodec2/unittest/tlspsens Executable file
View File

Binary file not shown.

127
libs/libcodec2/unittest/tlspsens.c Normal file
View File

@ -0,0 +1,127 @@
/*---------------------------------------------------------------------------*\
FILE........: tlspsens.c
AUTHOR......: David Rowe
DATE CREATED: 31 May 2012
Testing bit error sensitivity of LSP bits, first step in devising an unequal
error protection scheme.
\*---------------------------------------------------------------------------*/
/*
Copyright (C) 2012 David Rowe
All rights reserved.
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License version 2, as
published by the Free Software Foundation. This program is
distributed in the hope that it will be useful, but WITHOUT ANY
WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
License for more details.
You should have received a copy of the GNU Lesser General Public License
along with this program; if not, see <http://www.gnu.org/licenses/>.
*/
#include <assert.h>
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
#include "defines.h"
#include "comp.h"
#include "codec2.h"
#include "quantise.h"
#include "interp.h"
#include "codec2_internal.h"
float run_a_test(char raw_file_name[], int bit_to_corrupt)
{
FILE *fin;
short buf[N];
struct CODEC2 *c2;
kiss_fft_cfg fft_fwd_cfg;
MODEL model;
float ak[LPC_ORD+1];
float lsps[LPC_ORD], e;
int lsp_indexes[LPC_ORD], found_bit;
float snr, snr_sum;
int frames, i, mask, index;
c2 = codec2_create(CODEC2_MODE_2400);
fft_fwd_cfg = kiss_fft_alloc(FFT_ENC, 0, NULL, NULL);
fin = fopen(raw_file_name, "rb");
assert(fin != NULL);
/* find bit we are corrupting */
found_bit = 0;
for(i=0; i<LSP_SCALAR_INDEXES; i++) {
if (!found_bit) {
if (bit_to_corrupt > lsp_bits(i))
bit_to_corrupt -= lsp_bits(i);
else {
index = i;
mask = (1 << bit_to_corrupt);
printf(" index: %d bit: %d mask: 0x%x ", index, bit_to_corrupt, mask);
found_bit = 1;
}
}
}
assert(found_bit == 1);
/* OK test a sample file, flipping bit */
snr_sum = 0.0;
frames = 0;
while(fread(buf, sizeof(short), N, fin) == N) {
analyse_one_frame(c2, &model, buf);
e = speech_to_uq_lsps(lsps, ak, c2->Sn, c2->w, LPC_ORD);
encode_lsps_scalar(lsp_indexes, lsps, LPC_ORD);
/* find and flip bit we are testing */
lsp_indexes[index] ^= mask;
/* decode LSPs and measure SNR */
decode_lsps_scalar(lsps, lsp_indexes, LPC_ORD);
check_lsp_order(lsps, LPC_ORD);
bw_expand_lsps(lsps, LPC_ORD);
lsp_to_lpc(lsps, ak, LPC_ORD);
aks_to_M2(fft_fwd_cfg, ak, LPC_ORD, &model, e, &snr, 0, 0, 1, 1, LPCPF_BETA, LPCPF_GAMMA);
snr_sum += snr;
frames++;
}
codec2_destroy(c2);
fclose(fin);
return snr_sum/frames;
}
int main(int argc, char *argv[]) {
int i;
int total_lsp_bits = 0;
float snr;
if (argc != 2) {
printf("usage: %s RawFile\n", argv[0]);
exit(1);
}
for(i=0; i<LPC_ORD; i++)
total_lsp_bits += lsp_bits(i);
for(i=0; i<total_lsp_bits; i++) {
snr = run_a_test(argv[1], i);
printf("%d %5.2f\n", i, snr);
}
return 0;
}

BIN
libs/libcodec2/unittest/tprede Executable file
View File

Binary file not shown.

53
libs/libcodec2/unittest/tprede.c Normal file
View File

@ -0,0 +1,53 @@
/*
tpre_de.c
David Rowe
Sep 24 2012
Unit test to generate the combined impulse response of pre & de-emphasis filters.
pl("../unittest/out48.raw",1,3000)
pl("../unittest/out8.raw",1,3000)
Listening to it also shows up anything nasty:
$ play -s -2 -r 48000 out48.raw
$ play -s -2 -r 8000 out8.raw
*/
#include <assert.h>
#include <math.h>
#include <stdlib.h>
#include <stdio.h>
#include "lpc.h"
#define N 10
#define F 10
int main() {
FILE *fprede;
float Sn[N], Sn_pre[N], Sn_de[N];
float pre_mem = 0.0, de_mem = 0.0;
int i, f;
fprede = fopen("prede.txt", "wt");
assert(fprede != NULL);
for(i=0; i<N; i++)
Sn[i] = 0.0;
Sn[0]= 1.0;
for(f=0; f<F; f++) {
pre_emp(Sn_pre, Sn, &pre_mem, N);
de_emp(Sn_de, Sn_pre, &de_mem, N);
for(i=0; i<N; i++) {
fprintf(fprede, "%f\n", Sn_de[i]);
}
Sn[0] = 0.0;
}
fclose(fprede);
return 0;
}

64
libs/libcodec2/unittest/tsrc.c Normal file
View File

@ -0,0 +1,64 @@
/*
tsrc.c
David Rowe
Sat Nov 3 2012
Unit test for libresample code.
*/
#include <assert.h>
#include <math.h>
#include <stdlib.h>
#include <stdio.h>
#include <samplerate.h>
#define N8 160 /* processing buffer size at 8 kHz */
#define N48 ((int)N8*(48000/8000)) /* buf size assuming 48k max sample rate */
int main(int argc, char *argv[]) {
FILE *f8k, *fout;
short in8k_short[N8];
float in8k[N8];
float out[N48];
short out_short[N48];
SRC_STATE *src;
SRC_DATA data;
int error;
if (argc != 4) {
printf("usage %s inputRawFile OutputRawFile OutputSamplerate\n", argv[0]);
exit(0);
}
f8k = fopen(argv[1], "rb");
assert(f8k != NULL);
fout = fopen(argv[2], "wb");
assert(fout != NULL);
src = src_new(SRC_SINC_FASTEST, 1, &error);
assert(src != NULL);
data.data_in = in8k;
data.data_out = out;
data.input_frames = N8;
data.output_frames = N48;
data.end_of_input = 0;
data.src_ratio = atof(argv[3])/8000;
printf("%f\n", data.src_ratio);
while(fread(in8k_short, sizeof(short), N8, f8k) == N8) {
src_short_to_float_array(in8k_short, in8k, N8);
src_process(src, &data);
printf("%d %d\n", (int)data.output_frames , (int)data.output_frames_gen);
assert(data.output_frames_gen <= N48);
src_float_to_short_array(out, out_short, data.output_frames_gen);
fwrite(out_short, sizeof(short), data.output_frames_gen, fout);
}
fclose(fout);
fclose(f8k);
return 0;
}

BIN
libs/libcodec2/unittest/vq_train_jvm Executable file
View File

Binary file not shown.

486
libs/libcodec2/unittest/vq_train_jvm.c Executable file
View File

@ -0,0 +1,486 @@
/*---------------------------------------------------------------------------*\
FILE........: vq_train_jvm.c
AUTHOR......: Jean-Marc Valin
DATE CREATED: 21 Jan 2012
Multi-stage Vector Quantoser training program developed by Jean-Marc at
linux.conf.au 2012. Minor mods by David Rowe
\*---------------------------------------------------------------------------*/
/*
Copyright (C) 2012 Jean-Marc Valin
All rights reserved.
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License version 2, as
published by the Free Software Foundation. This program is
distributed in the hope that it will be useful, but WITHOUT ANY
WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
License for more details.
You should have received a copy of the GNU Lesser General Public License
along with this program; if not, see <http://www.gnu.org/licenses/>.
*/
#ifdef VALGRIND
#include <valgrind/memcheck.h>
#endif
#include <assert.h>
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#define MIN(a,b) ((a)<(b)?(a):(b))
#define COEF 0.0f
#define MAX_ENTRIES 16384
void compute_weights(const float *x, float *w, int ndim)
{
int i;
w[0] = MIN(x[0], x[1]-x[0]);
for (i=1;i<ndim-1;i++)
w[i] = MIN(x[i]-x[i-1], x[i+1]-x[i]);
w[ndim-1] = MIN(x[ndim-1]-x[ndim-2], M_PI-x[ndim-1]);
for (i=0;i<ndim;i++)
w[i] = 1./(.01+w[i]);
w[0]*=3;
w[1]*=2;
}
int find_nearest(const float *codebook, int nb_entries, float *x, int ndim, float *min_dist)
{
int i, j;
int nearest = 0;
*min_dist = 1E15;
for (i=0;i<nb_entries;i++)
{
float dist=0;
for (j=0;j<ndim;j++)
dist += (x[j]-codebook[i*ndim+j])*(x[j]-codebook[i*ndim+j]);
if (dist<*min_dist)
{
*min_dist = dist;
nearest = i;
}
}
return nearest;
}
int find_nearest_weighted(const float *codebook, int nb_entries, float *x, const float *w, int ndim)
{
int i, j;
float min_dist = 1e15;
int nearest = 0;
for (i=0;i<nb_entries;i++)
{
float dist=0;
for (j=0;j<ndim;j++)
dist += w[j]*(x[j]-codebook[i*ndim+j])*(x[j]-codebook[i*ndim+j]);
if (dist<min_dist)
{
min_dist = dist;
nearest = i;
}
}
return nearest;
}
int quantize_lsp(const float *x, const float *codebook1, const float *codebook2,
const float *codebook3, int nb_entries, float *xq, int ndim)
{
int i, n1, n2, n3;
float err[ndim], err2[ndim], err3[ndim];
float w[ndim], w2[ndim], w3[ndim], min_dist;
w[0] = MIN(x[0], x[1]-x[0]);
for (i=1;i<ndim-1;i++)
w[i] = MIN(x[i]-x[i-1], x[i+1]-x[i]);
w[ndim-1] = MIN(x[ndim-1]-x[ndim-2], M_PI-x[ndim-1]);
/*
for (i=0;i<ndim;i++)
w[i] = 1./(.003+w[i]);
w[0]*=3;
w[1]*=2;*/
compute_weights(x, w, ndim);
for (i=0;i<ndim;i++)
err[i] = x[i]-COEF*xq[i];
n1 = find_nearest(codebook1, nb_entries, err, ndim, &min_dist);
for (i=0;i<ndim;i++)
{
xq[i] = COEF*xq[i] + codebook1[ndim*n1+i];
err[i] -= codebook1[ndim*n1+i];
}
for (i=0;i<ndim/2;i++)
{
err2[i] = err[2*i];
err3[i] = err[2*i+1];
w2[i] = w[2*i];
w3[i] = w[2*i+1];
}
n2 = find_nearest_weighted(codebook2, nb_entries, err2, w2, ndim/2);
n3 = find_nearest_weighted(codebook3, nb_entries, err3, w3, ndim/2);
for (i=0;i<ndim/2;i++)
{
xq[2*i] += codebook2[ndim*n2/2+i];
xq[2*i+1] += codebook3[ndim*n3/2+i];
}
return 0;
}
void split(float *codebook, int nb_entries, int ndim)
{
int i,j;
for (i=0;i<nb_entries;i++)
{
for (j=0;j<ndim;j++)
{
float delta = .01*(rand()/(float)RAND_MAX-.5);
codebook[i*ndim+j] += delta;
codebook[(i+nb_entries)*ndim+j] = codebook[i*ndim+j] - delta;
}
}
}
void update(float *data, int nb_vectors, float *codebook, int nb_entries, int ndim)
{
int i,j;
int count[nb_entries];
int nearest[nb_vectors];
float min_dist;
float total_min_dist = 0;
for (i=0;i<nb_entries;i++)
count[i] = 0;
for (i=0;i<nb_vectors;i++)
{
nearest[i] = find_nearest(codebook, nb_entries, data+i*ndim, ndim, &min_dist);
total_min_dist += min_dist;
}
for (i=0;i<nb_entries*ndim;i++)
codebook[i] = 0;
for (i=0;i<nb_vectors;i++)
{
int n = nearest[i];
count[n]++;
for (j=0;j<ndim;j++)
codebook[n*ndim+j] += data[i*ndim+j];
}
float w2=0;
for (i=0;i<nb_entries;i++)
{
for (j=0;j<ndim;j++)
codebook[i*ndim+j] *= (1./count[i]);
w2 += (count[i]/(float)nb_vectors)*(count[i]/(float)nb_vectors);
}
fprintf(stderr, "%f / %d var = %f\n", 1./w2, nb_entries, total_min_dist/nb_vectors );
}
void update_weighted(float *data, float *weight, int nb_vectors, float *codebook, int nb_entries, int ndim)
{
int i,j;
float count[MAX_ENTRIES][ndim];
int nearest[nb_vectors];
for (i=0;i<nb_entries;i++)
for (j=0;j<ndim;j++)
count[i][j] = 0;
for (i=0;i<nb_vectors;i++)
{
nearest[i] = find_nearest_weighted(codebook, nb_entries, data+i*ndim, weight+i*ndim, ndim);
}
for (i=0;i<nb_entries*ndim;i++)
codebook[i] = 0;
for (i=0;i<nb_vectors;i++)
{
int n = nearest[i];
for (j=0;j<ndim;j++)
{
float w = sqrt(weight[i*ndim+j]);
count[n][j]+=w;
codebook[n*ndim+j] += w*data[i*ndim+j];
}
}
//float w2=0;
for (i=0;i<nb_entries;i++)
{
for (j=0;j<ndim;j++)
codebook[i*ndim+j] *= (1./count[i][j]);
//w2 += (count[i]/(float)nb_vectors)*(count[i]/(float)nb_vectors);
}
//fprintf(stderr, "%f / %d\n", 1./w2, nb_entries);
}
void vq_train(float *data, int nb_vectors, float *codebook, int nb_entries, int ndim)
{
int i, j, e;
e = 1;
for (j=0;j<ndim;j++)
codebook[j] = 0;
for (i=0;i<nb_vectors;i++)
for (j=0;j<ndim;j++)
codebook[j] += data[i*ndim+j];
for (j=0;j<ndim;j++)
codebook[j] *= (1./nb_vectors);
while (e< nb_entries)
{
split(codebook, e, ndim);
fprintf(stderr, "%d\n", e);
e<<=1;
for (j=0;j<ndim;j++)
update(data, nb_vectors, codebook, e, ndim);
}
}
void vq_train_weighted(float *data, float *weight, int nb_vectors, float *codebook, int nb_entries, int ndim)
{
int i, j, e;
e = 1;
for (j=0;j<ndim;j++)
codebook[j] = 0;
for (i=0;i<nb_vectors;i++)
for (j=0;j<ndim;j++)
codebook[j] += data[i*ndim+j];
for (j=0;j<ndim;j++)
codebook[j] *= (1./nb_vectors);
while (e<nb_entries)
{
split(codebook, e, ndim);
fprintf(stderr, "%d\n", e);
e<<=1;
for (j=0;j<ndim;j++)
update_weighted(data, weight, nb_vectors, codebook, e, ndim);
}
}
int main(int argc, char **argv)
{
int i,j;
FILE *ftrain;
int nb_vectors, nb_entries, ndim;
float *data, *pred, *codebook, *codebook2, *codebook3;
float *weight, *weight2, *weight3;
float *delta, *delta2;
float tmp, err, min_dist, total_min_dist;
int ret;
char filename[256];
FILE *fcb;
printf("Jean-Marc Valin's Split VQ training program....\n");
if (argc != 5) {
printf("usage: %s TrainTextFile K(dimension) M(codebook size) VQFilesPrefix\n", argv[0]);
exit(1);
}
ndim = atoi(argv[2]);
nb_vectors = atoi(argv[3]);
nb_entries = atoi(argv[3]);
/* determine size of training file */
ftrain = fopen(argv[1],"rt"); assert(ftrain != NULL);
nb_vectors = 0;
while (1) {
if (feof(ftrain))
break;
for (j=0;j<ndim;j++)
{
ret = fscanf(ftrain, "%f ", &tmp);
}
nb_vectors++;
if ((nb_vectors % 1000) == 0)
printf("\r%d lines",nb_vectors);
}
rewind(ftrain);
printf("\nndim %d nb_vectors %d nb_entries %d\n", ndim, nb_vectors, nb_entries);
data = malloc(nb_vectors*ndim*sizeof(*data));
weight = malloc(nb_vectors*ndim*sizeof(*weight));
weight2 = malloc(nb_vectors*ndim*sizeof(*weight2));
weight3 = malloc(nb_vectors*ndim*sizeof(*weight3));
pred = malloc(nb_vectors*ndim*sizeof(*pred));
codebook = malloc(nb_entries*ndim*sizeof(*codebook));
codebook2 = malloc(nb_entries*ndim*sizeof(*codebook2));
codebook3 = malloc(nb_entries*ndim*sizeof(*codebook3));
for (i=0;i<nb_vectors;i++)
{
if (feof(ftrain))
break;
for (j=0;j<ndim;j++)
{
ret = fscanf(ftrain, "%f ", &data[i*ndim+j]);
}
}
nb_vectors = i;
#ifdef VALGRIND
VALGRIND_CHECK_MEM_IS_DEFINED(data, nb_entries*ndim);
#endif
/* determine weights for each training vector */
for (i=0;i<nb_vectors;i++)
{
compute_weights(data+i*ndim, weight+i*ndim, ndim);
for (j=0;j<ndim/2;j++)
{
weight2[i*ndim/2+j] = weight[i*ndim+2*j];
weight3[i*ndim/2+j] = weight[i*ndim+2*j+1];
}
}
/* 20ms (two frame gaps) initial predictor state */
for (i=0;i<ndim;i++) {
pred[i+ndim] = pred[i] = data[i] - M_PI*(i+1)/(ndim+1);
}
/* generate predicted data for training */
for (i=2;i<nb_vectors;i++)
{
for (j=0;j<ndim;j++)
pred[i*ndim+j] = data[i*ndim+j] - COEF*data[(i-2)*ndim+j];
}
#ifdef VALGRIND
VALGRIND_CHECK_MEM_IS_DEFINED(pred, nb_entries*ndim);
#endif
/* train first stage */
vq_train(pred, nb_vectors, codebook, nb_entries, ndim);
delta = malloc(nb_vectors*ndim*sizeof(*data));
err = 0;
total_min_dist = 0;
for (i=0;i<nb_vectors;i++)
{
int nearest = find_nearest(codebook, nb_entries, &pred[i*ndim], ndim, &min_dist);
total_min_dist += min_dist;
for (j=0;j<ndim;j++)
{
//delta[i*ndim+j] = data[i*ndim+j] - codebook[nearest*ndim+j];
//printf("%f ", delta[i*ndim+j]);
//err += (delta[i*ndim+j])*(delta[i*ndim+j]);
delta[i*ndim/2+j/2+(j&1)*nb_vectors*ndim/2] = pred[i*ndim+j] - codebook[nearest*ndim+j];
//printf("%f ", delta[i*ndim/2+j/2+(j&1)*nb_vectors*ndim/2]);
err += (delta[i*ndim/2+j/2+(j&1)*nb_vectors*ndim/2])*(delta[i*ndim/2+j/2+(j&1)*nb_vectors*ndim/2]);
}
//printf("\n");
}
fprintf(stderr, "Stage 1 LSP RMS error: %f\n", sqrt(err/nb_vectors/ndim));
fprintf(stderr, "Stage 1 LSP variance.: %f\n", total_min_dist/nb_vectors);
#if 1
vq_train(delta, nb_vectors, codebook2, nb_entries, ndim/2);
vq_train(delta+ndim*nb_vectors/2, nb_vectors, codebook3, nb_entries, ndim/2);
#else
vq_train_weighted(delta, weight2, nb_vectors, codebook2, nb_entries, ndim/2);
vq_train_weighted(delta+ndim*nb_vectors/2, weight3, nb_vectors, codebook3, nb_entries, ndim/2);
#endif
err = 0;
total_min_dist = 0;
delta2 = delta + nb_vectors*ndim/2;
for (i=0;i<nb_vectors;i++)
{
int n1, n2;
n1 = find_nearest(codebook2, nb_entries, &delta[i*ndim/2], ndim/2, &min_dist);
for (j=0;j<ndim/2;j++)
{
delta[i*ndim/2+j] = delta[i*ndim/2+j] - codebook2[n1*ndim/2+j];
err += (delta[i*ndim/2+j])*(delta[i*ndim/2+j]);
}
total_min_dist += min_dist;
n2 = find_nearest(codebook3, nb_entries, &delta2[i*ndim/2], ndim/2, &min_dist);
for (j=0;j<ndim/2;j++)
{
delta[i*ndim/2+j] = delta[i*ndim/2+j] - codebook2[n2*ndim/2+j];
err += (delta2[i*ndim/2+j])*(delta2[i*ndim/2+j]);
}
total_min_dist += min_dist;
}
fprintf(stderr, "Stage 2 LSP RMS error: %f\n", sqrt(err/nb_vectors/ndim));
fprintf(stderr, "Stage 2 LSP Variance.: %f\n", total_min_dist/nb_vectors);
float xq[ndim];
for (i=0;i<ndim;i++)
xq[i] = M_PI*(i+1)/(ndim+1);
for (i=0;i<nb_vectors;i++)
{
quantize_lsp(data+i*ndim, codebook, codebook2,
codebook3, nb_entries, xq, ndim);
/*for (j=0;j<ndim;j++)
printf("%f ", xq[j]);
printf("\n");*/
}
/* save output tables to text files */
sprintf(filename, "%s1.txt", argv[4]);
fcb = fopen(filename, "wt"); assert(fcb != NULL);
fprintf(fcb, "%d %d\n", ndim, nb_entries);
for (i=0;i<nb_entries;i++)
{
for (j=0;j<ndim;j++)
fprintf(fcb, "%f ", codebook[i*ndim+j]);
fprintf(fcb, "\n");
}
fclose(fcb);
sprintf(filename, "%s2.txt", argv[4]);
fcb = fopen(filename, "wt"); assert(fcb != NULL);
fprintf(fcb, "%d %d\n", ndim/2, nb_entries);
for (i=0;i<nb_entries;i++)
{
for (j=0;j<ndim/2;j++)
fprintf(fcb, "%f ", codebook2[i*ndim/2+j]);
fprintf(fcb, "\n");
}
fclose(fcb);
sprintf(filename, "%s3.txt", argv[4]);
fcb = fopen(filename, "wt"); assert(fcb != NULL);
fprintf(fcb, "%d %d\n", ndim/2, nb_entries);
for (i=0;i<nb_entries;i++)
{
for (j=0;j<ndim/2;j++)
fprintf(fcb, "%f ", codebook3[i*ndim/2+j]);
fprintf(fcb, "\n");
}
fclose(fcb);
return 0;
}

BIN
libs/libcodec2/unittest/vqtrainjnd Executable file
View File

Binary file not shown.

BIN
libs/libcodec2/unittest/vqtrainph Executable file
View File

Binary file not shown.

419
libs/libcodec2/unittest/vqtrainph.c Normal file
View File

@ -0,0 +1,419 @@
/*--------------------------------------------------------------------------*\
FILE........: vqtrainph.c
AUTHOR......: David Rowe
DATE CREATED: 27 July 2012
This program trains phase vector quantisers. Modified from
vqtrain.c
\*--------------------------------------------------------------------------*/
/*
Copyright (C) 2012 David Rowe
All rights reserved.
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License version 2, as
published by the Free Software Foundation. This program is
distributed in the hope that it will be useful, but WITHOUT ANY
WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
License for more details.
You should have received a copy of the GNU Lesser General Public License
along with this program; if not, see <http://www.gnu.org/licenses/>.
*/
/*-----------------------------------------------------------------------*\
INCLUDES
\*-----------------------------------------------------------------------*/
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <ctype.h>
#include <assert.h>
typedef struct {
float real;
float imag;
} COMP;
/*-----------------------------------------------------------------------* \
DEFINES
\*-----------------------------------------------------------------------*/
#define DELTAQ 0.01 /* quiting distortion */
#define MAX_STR 80 /* maximum string length */
#define PI 3.141592654
/*-----------------------------------------------------------------------*\
FUNCTION PROTOTYPES
\*-----------------------------------------------------------------------*/
void zero(COMP v[], int d);
void acc(COMP v1[], COMP v2[], int d);
void norm(COMP v[], int k);
int quantise(COMP cb[], COMP vec[], int d, int e, float *se);
void print_vec(COMP cb[], int d, int e);
/*-----------------------------------------------------------------------* \
MAIN
\*-----------------------------------------------------------------------*/
int main(int argc, char *argv[]) {
int d,e; /* dimension and codebook size */
COMP *vec; /* current vector */
COMP *cb; /* vector codebook */
COMP *cent; /* centroids for each codebook entry */
int *n; /* number of vectors in this interval */
int J; /* number of vectors in training set */
int ind; /* index of current vector */
float se; /* total squared error for this iteration */
float var; /* variance */
float var_1; /* previous variance */
float delta; /* improvement in distortion */
FILE *ftrain; /* file containing training set */
FILE *fvq; /* file containing vector quantiser */
int ret;
int i,j, finished, iterations;
float b; /* equivalent number of bits */
float improvement;
float sd_vec, sd_element, sd_theory, bits_theory;
int var_n;
/* Interpret command line arguments */
if (argc != 5) {
printf("usage: %s TrainFile D(dimension) E(number of entries) VQFile\n", argv[0]);
exit(1);
}
/* Open training file */
ftrain = fopen(argv[1],"rb");
if (ftrain == NULL) {
printf("Error opening training database file: %s\n",argv[1]);
exit(1);
}
/* determine k and m, and allocate arrays */
d = atoi(argv[2]);
e = atoi(argv[3]);
printf("\n");
printf("dimension D=%d number of entries E=%d\n", d, e);
vec = (COMP*)malloc(sizeof(COMP)*d);
cb = (COMP*)malloc(sizeof(COMP)*d*e);
cent = (COMP*)malloc(sizeof(COMP)*d*e);
n = (int*)malloc(sizeof(int)*e);
if (cb == NULL || cb == NULL || cent == NULL || vec == NULL) {
printf("Error in malloc.\n");
exit(1);
}
/* determine size of training set */
J = 0;
var_n = 0;
while(fread(vec, sizeof(COMP), d, ftrain) == (size_t)d) {
for(j=0; j<d; j++)
if ((vec[j].real != 0.0) && (vec[j].imag != 0.0))
var_n++;
J++;
}
printf("J=%d sparse vectors in training set, %d non-zero phases\n", J, var_n);
/* set up initial codebook state from samples of training set */
rewind(ftrain);
ret = fread(cb, sizeof(COMP), d*e, ftrain);
/* codebook can't have any zero phase angle entries, these need to be set to
zero angle so cmult used to find phase angle differences works */
for(i=0; i<d*e; i++)
if ((cb[i].real == 0.0) && (cb[i].imag == 0.0)) {
cb[i].real = 1.0;
cb[i].imag = 0.0;
}
//print_vec(cb, d, 1);
/* main loop */
printf("\n");
printf("Iteration delta var std dev\n");
printf("--------------------------------\n");
b = log10((float)e)/log10(2.0);
sd_theory = (PI/sqrt(3.0))*pow(2.0, -b/(float)d);
iterations = 0;
finished = 0;
delta = 0;
var_1 = 0.0;
do {
/* zero centroids */
for(i=0; i<e; i++) {
zero(&cent[i*d], d);
n[i] = 0;
}
/* quantise training set */
se = 0.0;
rewind(ftrain);
for(i=0; i<J; i++) {
ret = fread(vec, sizeof(COMP), d, ftrain);
ind = quantise(cb, vec, d, e, &se);
//printf("%d ", ind);
n[ind]++;
acc(&cent[ind*d], vec, d);
}
/* work out stats */
var = se/var_n;
sd_vec = sqrt(var);
/* we need to know dimension of cb (which varies from vector to vector)
to calc bits_theory. Maybe measure and use average dimension....
*/
//bits_theory = d*log10(PI/(sd_element*sqrt(3.0)))/log10(2.0);
//improvement = bits_theory - b;
//print_vec(cent, d, 1);
//print_vec(cb, d, 1);
iterations++;
if (iterations > 1) {
if (var > 0.0) {
delta = (var_1 - var)/var;
}
else
delta = 0;
if (delta < DELTAQ)
finished = 1;
}
if (!finished) {
/* determine new codebook from centroids */
for(i=0; i<e; i++) {
norm(&cent[i*d], d);
memcpy(&cb[i*d], &cent[i*d], d*sizeof(COMP));
}
}
printf("%2d %4.3f %4.3f %4.3f \n",iterations, delta, var, sd_vec);
var_1 = var;
} while (!finished);
//print_vec(cb, d, 1);
/* save codebook to disk */
fvq = fopen(argv[4],"wt");
if (fvq == NULL) {
printf("Error opening VQ file: %s\n",argv[4]);
exit(1);
}
fprintf(fvq,"%d %d\n",d,e);
for(j=0; j<e; j++) {
for(i=0; i<d; i++)
fprintf(fvq,"% 4.3f ", atan2(cb[j*d+i].imag, cb[j*d+i].real));
fprintf(fvq,"\n");
}
fclose(fvq);
return 0;
}
/*-----------------------------------------------------------------------*\
FUNCTIONS
\*-----------------------------------------------------------------------*/
void print_vec(COMP cb[], int d, int e)
{
int i,j;
for(j=0; j<e; j++) {
for(i=0; i<d; i++)
printf("%f %f ", cb[j*d+i].real, cb[j*d+i].imag);
printf("\n");
}
}
static COMP cconj(COMP a)
{
COMP res;
res.real = a.real;
res.imag = -a.imag;
return res;
}
static COMP cmult(COMP a, COMP b)
{
COMP res;
res.real = a.real*b.real - a.imag*b.imag;
res.imag = a.real*b.imag + a.imag*b.real;
return res;
}
static COMP cadd(COMP a, COMP b)
{
COMP res;
res.real = a.real + b.real;
res.imag = a.imag + b.imag;
return res;
}
/*---------------------------------------------------------------------------*\
FUNCTION....: zero()
AUTHOR......: David Rowe
DATE CREATED: 23/2/95
Zeros a vector of length d.
\*---------------------------------------------------------------------------*/
void zero(COMP v[], int d)
{
int i;
for(i=0; i<d; i++) {
v[i].real = 0.0;
v[i].imag = 0.0;
}
}
/*---------------------------------------------------------------------------*\
FUNCTION....: acc()
AUTHOR......: David Rowe
DATE CREATED: 23/2/95
Adds d dimensional vectors v1 to v2 and stores the result back
in v1. We add them like vectors on the complex plane, summing
the real and imag terms.
An unused entry in a sparse vector has both the real and imag
parts set to zero so won't affect the accumulation process.
\*---------------------------------------------------------------------------*/
void acc(COMP v1[], COMP v2[], int d)
{
int i;
for(i=0; i<d; i++)
v1[i] = cadd(v1[i], v2[i]);
}
/*---------------------------------------------------------------------------*\
FUNCTION....: norm()
AUTHOR......: David Rowe
DATE CREATED: 23/2/95
Normalises each element in d dimensional vector.
\*---------------------------------------------------------------------------*/
void norm(COMP v[], int d)
{
int i;
float mag;
for(i=0; i<d; i++) {
mag = sqrt(v[i].real*v[i].real + v[i].imag*v[i].imag);
if (mag == 0.0) {
/* can't have zero cb entries as cmult will break in quantise().
We effectively set sparese phases to an angle of 0. */
v[i].real = 1.0;
v[i].imag = 0.0;
}
else {
v[i].real /= mag;
v[i].imag /= mag;
}
}
}
/*---------------------------------------------------------------------------*\
FUNCTION....: quantise()
AUTHOR......: David Rowe
DATE CREATED: 23/2/95
Quantises vec by choosing the nearest vector in codebook cb, and
returns the vector index. The squared error of the quantised vector
is added to se.
Unused entries in sparse vectors are ignored.
\*---------------------------------------------------------------------------*/
int quantise(COMP cb[], COMP vec[], int d, int e, float *se)
{
float error; /* current error */
int besti; /* best index so far */
float best_error; /* best error so far */
int i,j;
int ignore;
COMP diff;
besti = 0;
best_error = 1E32;
for(j=0; j<e; j++) {
error = 0.0;
for(i=0; i<d; i++) {
ignore = (vec[i].real == 0.0) && (vec[i].imag == 0.0);
if (!ignore) {
diff = cmult(cb[j*d+i], cconj(vec[i]));
error += pow(atan2(diff.imag, diff.real), 2.0);
}
}
if (error < best_error) {
best_error = error;
besti = j;
}
}
*se += best_error;
return(besti);
}

Some files were not shown because too many files have changed in this diff Show More