ref: b64889efb02aab61b670ed228e5790c865b053e2
parent: b16609338307966dd259a871d3cc94a0294b51e0
author: Erik de Castro Lopo <erikd@miles>
date: Sun Oct 31 17:30:24 EST 2004
Sketch out skeleton of program fir_interp.
--- /dev/null
+++ b/FIR-POLY/Makefile
@@ -1,0 +1,40 @@
+TARGETS = fir_interp
+
+CXX = g++
+CXXFLAGS = -g -O3 -Wall -W -Werror -Wstrict-prototypes -Wmissing-prototypes
+
+GMATH_INC = $(shell pkg-config --cflags GMath)
+GMATH_LIB = $(shell pkg-config --libs GMath)
+
+FFTW_INC = $(shell pkg-config --cflags fftw3)
+FFTW_LIB = $(shell pkg-config --libs fftw3)
+
+
+
+all: $(TARGETS)
+
+clean:
+ rm -f $(TARGETS) *.o
+
+check: $(TARGETS)
+ ./fir_interp
+
+fir_interp: fir_interp.o mag_spectrum.o demin.o
+ $(CXX) $(CXXFLAGS) $+ $(GMATH_LIB) $(FFTW_LIB) -o $@
+
+#---------------------------------------------------------------------
+
+fir_interp.o : fir_interp.cc mag_spectrum.hh
+ $(CXX) $(CXXFLAGS) $(GMATH_INC) -c fir_interp.cc -o $@
+
+mag_spectrum.o : mag_spectrum.cc mag_spectrum.hh
+ $(CXX) $(CXXFLAGS) $(FFTW_INC) -c mag_spectrum.cc -o $@
+
+demin.o : demin.cc
+ $(CXX) $(CXXFLAGS) $(GMATH_INC) -c $+ -o $@
+
+# Do not edit or modify anything in this comment block.
+# The arch-tag line is a file identity tag for the GNU Arch
+# revision control system.
+#
+# arch-tag: 34b6b183-dca1-4085-9089-c8b7d5b09099
--- /dev/null
+++ b/FIR-POLY/demin.cc
@@ -1,0 +1,268 @@
+/*
+** Copyright (C) 1998-2004 Erik de Castro Lopo <erikd@mega-nerd.com>
+**
+** The copyright above and this notice must be preserved in all
+** copies of this source code. The copyright above does not
+** evidence any actual or intended publication of this source code.
+**
+** This is unpublished proprietary trade secret of Erik de Castro
+** Lopo. This source code may not be copied, disclosed,
+** distributed, demonstrated or licensed except as authorized by
+** Erik de Castro Lopo.
+**
+** No part of this program or publication may be reproduced,
+** transmitted, transcribed, stored in a retrieval system,
+** or translated into any language or computer language in any
+** form or by any means, electronic, mechanical, magnetic,
+** optical, chemical, manual, or otherwise, without the prior
+** written permission of Erik de Castro Lopo.
+*/
+
+#include <stdio.h>
+#include <unistd.h>
+
+#include <Minimize.hh>
+#include <PsuedoRand.hh>
+
+class DEmin
+{ private :+ u_int PopSize, FuncCalcs, Generations ;
+
+ PsuedoRand Rand ;
+
+ GMatrix *MainVector, *Vector2 ;
+ double *MainCost, *Cost2 ;
+ double *ScaleFactor, *ScaleFactor2 ;
+ double *Crossover, *Crossover2 ;
+
+ MinFunc Function ;
+
+ public :
+ DEmin (void) ;
+ ~DEmin (void) ;
+
+ void Setup (MinFunc func, const GMatrix& start, u_int popsize, double range, u_long seed) ;
+
+ GMatrix Minimize (double tol, u_int maxgen) ;
+
+ void PrintStats (void) ;
+
+} ; // class DEmin
+
+//======================================================================================
+
+inline
+void Swap (GMatrix **a, GMatrix **b)
+{ GMatrix *temp ;+ temp = *a ; *a = *b ; *b = temp ;
+} ; //
+
+inline void
+Swap (double **a, double **b)
+{ double *temp ;+ temp = *a ; *a = *b ; *b = temp ;
+} ; // Swap
+
+inline double
+Max (double a, double b)
+{ return (a > b) ? a : b ;+} ; // Max
+
+inline double
+Min (double a, double b)
+{ return (a < b) ? a : b ;+} ; // Max
+
+//======================================================================================
+
+double
+MinDiffEvol (MinFunc Function, GMatrix& start, double Ftol, double range, u_long seedval)
+{ DEmin *pfmin ;+ u_int dimensions ;
+
+ dimensions = start.GetElements () ;
+
+ pfmin = new DEmin ;
+
+ pfmin -> Setup (Function, start, dimensions * 5, range, seedval) ;
+
+ Ftol = (Ftol < 1e-8) ? 1e-8 : Ftol ;
+
+ start = pfmin -> Minimize (Ftol, dimensions * 5000) ;
+
+ delete pfmin ;
+
+ return Function (start) ;
+} ; // MinDiffEvol
+
+
+//======================================================================================
+
+DEmin::DEmin (void)
+{ PopSize = 0 ;+ FuncCalcs = 0 ;
+ Generations = 0 ;
+ MainVector = NULL ;
+ MainCost = NULL ;
+ ScaleFactor = NULL ;
+ Crossover = NULL ;
+ Vector2 = NULL ;
+ Cost2 = NULL ;
+ ScaleFactor2 = NULL ;
+ Crossover2 = NULL ;
+} ; // DEmin constructor
+
+DEmin::~DEmin (void)
+{ if (PopSize)+ { delete [] MainVector ;+ delete MainCost ;
+ delete ScaleFactor ;
+ delete Crossover ;
+ delete [] Vector2 ;
+ delete Cost2 ;
+ delete ScaleFactor2 ;
+ delete Crossover2 ;
+ } ;
+ PopSize = 0 ;
+} ; // DEmin destructor
+
+void
+DEmin::Setup (MinFunc func, const GMatrix& start, u_int popsize, double range, u_long seed)
+{ u_int k ;+ u_short rows, cols ;
+
+ Generations = 0 ;
+ Function = func ;
+
+ PopSize = (popsize < 2) ? popsize + 2 : popsize ;
+
+ MainVector = new GMatrix [PopSize] ;
+ MainCost = new double [PopSize] ;
+ ScaleFactor = new double [PopSize] ;
+ Crossover = new double [PopSize] ;
+ Vector2 = new GMatrix [PopSize] ;
+ Cost2 = new double [PopSize] ;
+ ScaleFactor2 = new double [PopSize] ;
+ Crossover2 = new double [PopSize] ;
+
+ Rand.Seed (seed + 1234567) ;
+
+ rows = start.GetRows () ;
+ cols = start.GetCols () ;
+
+ GMatrix::PsuedoRandSeed (123123) ;
+
+ if (IsComplex (start))
+ for (k = 0 ; k < PopSize ; k++)
+ { MainVector [k].SetSizeComplex (rows, cols) ;+ MainVector [k].PsuedoRandomize () ;
+ MainVector [k] = (2.0 * range) * (MainVector [k] - Complex (0.5, 0.5)) ;
+ }
+ else
+ for (k = 0 ; k < PopSize ; k++)
+ { MainVector [k].SetSizeDouble (rows, cols) ;+ MainVector [k].PsuedoRandomize () ;
+ MainVector [k] = (2.0 * range) * (MainVector [k] - 0.5) ;
+ } ;
+
+ MainVector [0] = start ; // Must include given vector.
+
+ for (k = 0 ; k < PopSize ; k++)
+ { MainCost [k] = Function (MainVector [k]) ;+ FuncCalcs ++ ;
+ ScaleFactor [k] = 1.2 * Rand.UDDouble () ;
+ Crossover [k] = Rand.UDDouble () ;
+ } ;
+
+ return ;
+} ; // DEmin::Setup
+
+GMatrix
+DEmin::Minimize (double Ftol, u_int maxgen)
+{ u_int k, a, b, c ;+ GMatrix temp ;
+ double scalefactor, crossover, cost, improvement, best ;
+
+ improvement = 100 * Ftol ;
+ best = 10e10 ;
+
+ while (Generations < maxgen && best > Ftol)
+ { Swap (&MainVector, &Vector2) ;+ Swap (&MainCost, &Cost2) ;
+ Swap (&ScaleFactor, &ScaleFactor2) ;
+ Swap (&Crossover, &Crossover2) ;
+
+ improvement = 0.0 ;
+
+ for (k = 0 ; k < PopSize ; k++)
+ { do+ a = u_int (double (PopSize) * Rand.UDDouble ()) ;
+ while (a == k) ;
+ do
+ b = u_int (double (PopSize) * Rand.UDDouble ()) ;
+ while (b == k || b == a) ;
+ do
+ c = u_int (double (PopSize) * Rand.UDDouble ()) ;
+ while (c == k || c == a || c == b) ;
+
+if (a >= PopSize) { printf ("a - very bad.\n") ; exit (0) ; } ;+if (b >= PopSize) { printf ("b - very bad.\n") ; exit (0) ; } ;+if (c >= PopSize) { printf ("c - very bad.\n") ; exit (0) ; } ;+
+ scalefactor = 0.5 * (ScaleFactor2 [b] + ScaleFactor2 [c]) ;
+
+ temp = Vector2 [a] + scalefactor * (Vector2 [b] - Vector2 [c]) ;
+
+ crossover = 0.5 * (Crossover2 [k] + Crossover2 [a]) ;
+
+ temp.CrossOver (Vector2 [k], crossover) ;
+
+ cost = Function (temp) ;
+ FuncCalcs ++ ;
+
+ if (cost < Cost2 [k])
+ { MainVector [k] = temp ;+ MainCost [k] = cost ;
+ ScaleFactor [k] = scalefactor ;
+ Crossover [k] = crossover ;
+ improvement = 0.5 * (improvement + Max (improvement, Cost2 [k] - cost)) ;
+ best = Min (best, cost) ;
+ }
+ else
+ { MainVector [k] = Vector2 [k] ;+ MainCost [k] = Cost2 [k] ;
+ ScaleFactor [k] = ScaleFactor2 [k] ;
+ Crossover [k] = Crossover2 [k] ;
+ } ; // if
+ } ; // for k
+
+ Generations ++ ;
+ } ; // while
+
+ // Now find best vector.
+
+ for (a = 0, k = 1 ; k < PopSize ; k++)
+ if (MainCost [k] < MainCost [a])
+ a = k ;
+ temp = MainVector [a] ;
+
+ return temp ;
+} ; // DEmin::Minimize
+
+void
+DEmin::PrintStats (void)
+{ printf ("Generations : %8d\n", Generations) ;+ printf ("Function evaluations : %8d\n", FuncCalcs) ;+ putchar ('\n') ;+} ; // DEmin::PrintStats
+
+//======================================================================================
+
+
+
+// Do not edit or modify anything in this comment block.
+// The arch-tag line is a file identity tag for the GNU Arch
+// revision control system.
+//
+// arch-tag: e874ad11-7294-4757-9370-fd1e6e1261d0
+
--- /dev/null
+++ b/FIR-POLY/fir_interp.cc
@@ -1,0 +1,109 @@
+/*
+** Copyright (C) 2004 Erik de Castro Lopo <erikd@mega-nerd.com>
+**
+** The copyright above and this notice must be preserved in all
+** copies of this source code. The copyright above does not
+** evidence any actual or intended publication of this source code.
+**
+** This is unpublished proprietary trade secret of Erik de Castro
+** Lopo. This source code may not be copied, disclosed,
+** distributed, demonstrated or licensed except as authorized by
+** Erik de Castro Lopo.
+**
+** No part of this program or publication may be reproduced,
+** transmitted, transcribed, stored in a retrieval system,
+** or translated into any language or computer language in any
+** form or by any means, electronic, mechanical, magnetic,
+** optical, chemical, manual, or otherwise, without the prior
+** written permission of Erik de Castro Lopo.
+*/
+
+
+#include <cstdio>
+#include <cstdlib>
+#include <ctime>
+#include <cstring>
+#include <complex>
+
+#include <GMatrix.hh>
+#include <Minimize.hh>
+
+#include "mag_spectrum.hh"
+
+
+#define ARRAY_LEN(x) (sizeof (x) / sizeof ((x) [0]))
+
+/*
+** Number of half cycles of impulse response to the left and right
+** of the center.
+** left_half_offset <= right_half_cycles
+*/
+
+// static int left_half_cycles = 5 ;
+// static int right_half_cycles = 15 ;
+
+
+typedef struct
+{+ double sinc_time_fudge ;
+
+ /* Window coefficients (all in range [0, 1]). */
+ double left_ak [5], right_ak [5] ;
+
+} FIR_PARAMS ;
+
+static double fir_error (const GMatrix& gm) ;
+
+int
+main (void)
+{ static double data [32] ;+ GMatrix gm_start ;
+ double error ;
+ int k, param_count ;
+
+ param_count = sizeof (FIR_PARAMS) / sizeof (double) ;
+
+ srandom (time (NULL)) ;
+ for (k = 0 ; k < param_count ; k++)
+ data [k] = 3.0 * ((1.0 * random ()) / INT_MAX - 0.5) ;
+
+ gm_start = GMatrix (param_count, 1, data) ;
+
+ fir_error (gm_start) ;
+
+ do
+ error = MinDiffEvol (fir_error, gm_start, 1e-15, 2.0, random ()) ;
+ while (error > 10.0) ;
+
+ printf ("error : %f\n", error) ;+
+ return 0 ;
+} /* main */
+
+/*==============================================================================
+*/
+
+static double
+fir_error (const GMatrix& gm)
+{ FIR_PARAMS params ;+ double error = 0.0 ;
+ unsigned param_count ;
+
+
+ param_count = sizeof (params) / sizeof (double) ;
+
+ if (gm.GetData (param_count, (double*) ¶ms) != param_count)
+ { printf ("\n\nError : GetData should return 2.\n") ;+ exit (1) ;
+ } ;
+
+ return error ;
+} /* fir_error */
+
+/*
+** Do not edit or modify anything in this comment block.
+** The following line is a file identity tag for the GNU Arch
+** revision control system.
+**
+** arch-tag: 3ce7ca6f-394e-432c-aeaa-f228afc05afd
+*/
--- /dev/null
+++ b/FIR-POLY/mag_spectrum.cc
@@ -1,0 +1,62 @@
+/*
+** Copyright (C) 2004 Erik de Castro Lopo <erikd@mega-nerd.com>
+**
+** The copyright above and this notice must be preserved in all
+** copies of this source code. The copyright above does not
+** evidence any actual or intended publication of this source code.
+**
+** This is unpublished proprietary trade secret of Erik de Castro
+** Lopo. This source code may not be copied, disclosed,
+** distributed, demonstrated or licensed except as authorized by
+** Erik de Castro Lopo.
+**
+** No part of this program or publication may be reproduced,
+** transmitted, transcribed, stored in a retrieval system,
+** or translated into any language or computer language in any
+** form or by any means, electronic, mechanical, magnetic,
+** optical, chemical, manual, or otherwise, without the prior
+** written permission of Erik de Castro Lopo.
+*/
+
+#include <cstdio>
+#include <cstdlib>
+#include <cstring>
+#include <cmath>
+
+#include <fftw3.h>
+
+void
+mag_spectrum (double *input, int len, double *magnitude, double scale)
+{ static fftw_plan plan = NULL ;+
+ double maxval ;
+ int k ;
+
+ if (input == NULL || magnitude == NULL)
+ return ;
+
+ if (plan == NULL)
+ { plan = fftw_plan_r2r_1d (len, input, magnitude, FFTW_R2HC, FFTW_ESTIMATE | FFTW_PRESERVE_INPUT) ;+ if (plan == NULL)
+ { printf ("%s : line %d : create plan failed.\n", __FILE__, __LINE__) ;+ exit (1) ;
+ } ;
+ } ;
+
+ fftw_execute (plan) ;
+
+ /* (k < N/2 rounded up) */
+ maxval = 0.0 ;
+ for (k = 0 ; k < len / 2 ; k++)
+ magnitude [k] = scale * sqrt (magnitude [k] * magnitude [k] + magnitude [len - k - 1] * magnitude [len - k - 1]) ;
+
+ memset (magnitude + len / 2, 0, len / 2 * sizeof (magnitude [0])) ;
+
+ return ;
+} /* mag_spectrum */
+
+// Do not edit or modify anything in this comment block.
+// The following line is a file identity tag for the GNU Arch
+// revision control system.
+//
+// arch-tag: 58aa43ee-b8e1-4178-8861-7621399fa049
--- /dev/null
+++ b/FIR-POLY/mag_spectrum.hh
@@ -1,0 +1,27 @@
+/*
+** Copyright (C) 2004 Erik de Castro Lopo <erikd@mega-nerd.com>
+**
+** The copyright above and this notice must be preserved in all
+** copies of this source code. The copyright above does not
+** evidence any actual or intended publication of this source code.
+**
+** This is unpublished proprietary trade secret of Erik de Castro
+** Lopo. This source code may not be copied, disclosed,
+** distributed, demonstrated or licensed except as authorized by
+** Erik de Castro Lopo.
+**
+** No part of this program or publication may be reproduced,
+** transmitted, transcribed, stored in a retrieval system,
+** or translated into any language or computer language in any
+** form or by any means, electronic, mechanical, magnetic,
+** optical, chemical, manual, or otherwise, without the prior
+** written permission of Erik de Castro Lopo.
+*/
+
+void mag_spectrum (double *input, int len, double *magnitude, double scale) ;
+
+// Do not edit or modify anything in this comment block.
+// The following line is a file identity tag for the GNU Arch
+// revision control system.
+//
+// arch-tag: cc5db1c8-4c67-49f7-9c2f-6c5021a69970