Sketch out skeleton of program fir_interp.

Erik de Castro Lopo committed 22 years ago

Commit b64889e1 parent b1660935 files · +506 −0Browse files

Unified · Split

5 files

FIR-POLY/Makefile +40 −0
FIR-POLY/demin.cc +268 −0
FIR-POLY/fir_interp.cc +109 −0
FIR-POLY/mag_spectrum.cc +62 −0
FIR-POLY/mag_spectrum.hh +27 −0

@@ -0,0 +1,40 @@
	1	+TARGETS = fir_interp
	2	+
	3	+CXX = g++
	4	+CXXFLAGS = -g -O3 -Wall -W -Werror -Wstrict-prototypes -Wmissing-prototypes
	5	+
	6	+GMATH_INC = $(shell pkg-config --cflags GMath)
	7	+GMATH_LIB = $(shell pkg-config --libs GMath)
	8	+
	9	+FFTW_INC = $(shell pkg-config --cflags fftw3)
	10	+FFTW_LIB = $(shell pkg-config --libs fftw3)
	11	+
	12	+
	13	+
	14	+all: $(TARGETS)
	15	+
	16	+clean:
	17	+ rm -f $(TARGETS) *.o
	18	+
	19	+check: $(TARGETS)
	20	+ ./fir_interp
	21	+
	22	+fir_interp: fir_interp.o mag_spectrum.o demin.o
	23	+ $(CXX) $(CXXFLAGS) $+ $(GMATH_LIB) $(FFTW_LIB) -o $@
	24	+
	25	+#---------------------------------------------------------------------
	26	+
	27	+fir_interp.o : fir_interp.cc mag_spectrum.hh
	28	+ $(CXX) $(CXXFLAGS) $(GMATH_INC) -c fir_interp.cc -o $@
	29	+
	30	+mag_spectrum.o : mag_spectrum.cc mag_spectrum.hh
	31	+ $(CXX) $(CXXFLAGS) $(FFTW_INC) -c mag_spectrum.cc -o $@
	32	+
	33	+demin.o : demin.cc
	34	+ $(CXX) $(CXXFLAGS) $(GMATH_INC) -c $+ -o $@
	35	+
	36	+# Do not edit or modify anything in this comment block.
	37	+# The arch-tag line is a file identity tag for the GNU Arch
	38	+# revision control system.
	39	+#
	40	+# arch-tag: 34b6b183-dca1-4085-9089-c8b7d5b09099

FIR-POLY/demin.cc+268−0

@@ -0,0 +1,268 @@
	1	+/*
	2	+** Copyright (C) 1998-2004 Erik de Castro Lopo <erikd@mega-nerd.com>
	3	+**
	4	+** The copyright above and this notice must be preserved in all
	5	+** copies of this source code. The copyright above does not
	6	+** evidence any actual or intended publication of this source code.
	7	+**
	8	+** This is unpublished proprietary trade secret of Erik de Castro
	9	+** Lopo. This source code may not be copied, disclosed,
	10	+** distributed, demonstrated or licensed except as authorized by
	11	+** Erik de Castro Lopo.
	12	+**
	13	+** No part of this program or publication may be reproduced,
	14	+** transmitted, transcribed, stored in a retrieval system,
	15	+** or translated into any language or computer language in any
	16	+** form or by any means, electronic, mechanical, magnetic,
	17	+** optical, chemical, manual, or otherwise, without the prior
	18	+** written permission of Erik de Castro Lopo.
	19	+*/
	20	+
	21	+#include <stdio.h>
	22	+#include <unistd.h>
	23	+
	24	+#include <Minimize.hh>
	25	+#include <PsuedoRand.hh>
	26	+
	27	+class DEmin
	28	+{ private :
	29	+ u_int PopSize, FuncCalcs, Generations ;
	30	+
	31	+ PsuedoRand Rand ;
	32	+
	33	+ GMatrix MainVector, Vector2 ;
	34	+ double MainCost, Cost2 ;
	35	+ double ScaleFactor, ScaleFactor2 ;
	36	+ double Crossover, Crossover2 ;
	37	+
	38	+ MinFunc Function ;
	39	+
	40	+ public :
	41	+ DEmin (void) ;
	42	+ ~DEmin (void) ;
	43	+
	44	+ void Setup (MinFunc func, const GMatrix& start, u_int popsize, double range, u_long seed) ;
	45	+
	46	+ GMatrix Minimize (double tol, u_int maxgen) ;
	47	+
	48	+ void PrintStats (void) ;
	49	+
	50	+} ; // class DEmin
	51	+
	52	+//======================================================================================
	53	+
	54	+inline
	55	+void Swap (GMatrix a, GMatrix b)
	56	+{ GMatrix *temp ;
	57	+ temp = a ; a = b ; b = temp ;
	58	+} ; //
	59	+
	60	+inline void
	61	+Swap (double a, double b)
	62	+{ double *temp ;
	63	+ temp = a ; a = b ; b = temp ;
	64	+} ; // Swap
	65	+
	66	+inline double
	67	+Max (double a, double b)
	68	+{ return (a > b) ? a : b ;
	69	+} ; // Max
	70	+
	71	+inline double
	72	+Min (double a, double b)
	73	+{ return (a < b) ? a : b ;
	74	+} ; // Max
	75	+
	76	+//======================================================================================
	77	+
	78	+double
	79	+MinDiffEvol (MinFunc Function, GMatrix& start, double Ftol, double range, u_long seedval)
	80	+{ DEmin *pfmin ;
	81	+ u_int dimensions ;
	82	+
	83	+ dimensions = start.GetElements () ;
	84	+
	85	+ pfmin = new DEmin ;
	86	+
	87	+ pfmin -> Setup (Function, start, dimensions * 5, range, seedval) ;
	88	+
	89	+ Ftol = (Ftol < 1e-8) ? 1e-8 : Ftol ;
	90	+
	91	+ start = pfmin -> Minimize (Ftol, dimensions * 5000) ;
	92	+
	93	+ delete pfmin ;
	94	+
	95	+ return Function (start) ;
	96	+} ; // MinDiffEvol
	97	+
	98	+
	99	+//======================================================================================
	100	+
	101	+DEmin::DEmin (void)
	102	+{ PopSize = 0 ;
	103	+ FuncCalcs = 0 ;
	104	+ Generations = 0 ;
	105	+ MainVector = NULL ;
	106	+ MainCost = NULL ;
	107	+ ScaleFactor = NULL ;
	108	+ Crossover = NULL ;
	109	+ Vector2 = NULL ;
	110	+ Cost2 = NULL ;
	111	+ ScaleFactor2 = NULL ;
	112	+ Crossover2 = NULL ;
	113	+} ; // DEmin constructor
	114	+
	115	+DEmin::~DEmin (void)
	116	+{ if (PopSize)
	117	+ { delete [] MainVector ;
	118	+ delete MainCost ;
	119	+ delete ScaleFactor ;
	120	+ delete Crossover ;
	121	+ delete [] Vector2 ;
	122	+ delete Cost2 ;
	123	+ delete ScaleFactor2 ;
	124	+ delete Crossover2 ;
	125	+ } ;
	126	+ PopSize = 0 ;
	127	+} ; // DEmin destructor
	128	+
	129	+void
	130	+DEmin::Setup (MinFunc func, const GMatrix& start, u_int popsize, double range, u_long seed)
	131	+{ u_int k ;
	132	+ u_short rows, cols ;
	133	+
	134	+ Generations = 0 ;
	135	+ Function = func ;
	136	+
	137	+ PopSize = (popsize < 2) ? popsize + 2 : popsize ;
	138	+
	139	+ MainVector = new GMatrix [PopSize] ;
	140	+ MainCost = new double [PopSize] ;
	141	+ ScaleFactor = new double [PopSize] ;
	142	+ Crossover = new double [PopSize] ;
	143	+ Vector2 = new GMatrix [PopSize] ;
	144	+ Cost2 = new double [PopSize] ;
	145	+ ScaleFactor2 = new double [PopSize] ;
	146	+ Crossover2 = new double [PopSize] ;
	147	+
	148	+ Rand.Seed (seed + 1234567) ;
	149	+
	150	+ rows = start.GetRows () ;
	151	+ cols = start.GetCols () ;
	152	+
	153	+ GMatrix::PsuedoRandSeed (123123) ;
	154	+
	155	+ if (IsComplex (start))
	156	+ for (k = 0 ; k < PopSize ; k++)
	157	+ { MainVector [k].SetSizeComplex (rows, cols) ;
	158	+ MainVector [k].PsuedoRandomize () ;
	159	+ MainVector [k] = (2.0 * range) * (MainVector [k] - Complex (0.5, 0.5)) ;
	160	+ }
	161	+ else
	162	+ for (k = 0 ; k < PopSize ; k++)
	163	+ { MainVector [k].SetSizeDouble (rows, cols) ;
	164	+ MainVector [k].PsuedoRandomize () ;
	165	+ MainVector [k] = (2.0 * range) * (MainVector [k] - 0.5) ;
	166	+ } ;
	167	+
	168	+ MainVector [0] = start ; // Must include given vector.
	169	+
	170	+ for (k = 0 ; k < PopSize ; k++)
	171	+ { MainCost [k] = Function (MainVector [k]) ;
	172	+ FuncCalcs ++ ;
	173	+ ScaleFactor [k] = 1.2 * Rand.UDDouble () ;
	174	+ Crossover [k] = Rand.UDDouble () ;
	175	+ } ;
	176	+
	177	+ return ;
	178	+} ; // DEmin::Setup
	179	+
	180	+GMatrix
	181	+DEmin::Minimize (double Ftol, u_int maxgen)
	182	+{ u_int k, a, b, c ;
	183	+ GMatrix temp ;
	184	+ double scalefactor, crossover, cost, improvement, best ;
	185	+
	186	+ improvement = 100 * Ftol ;
	187	+ best = 10e10 ;
	188	+
	189	+ while (Generations < maxgen && best > Ftol)
	190	+ { Swap (&MainVector, &Vector2) ;
	191	+ Swap (&MainCost, &Cost2) ;
	192	+ Swap (&ScaleFactor, &ScaleFactor2) ;
	193	+ Swap (&Crossover, &Crossover2) ;
	194	+
	195	+ improvement = 0.0 ;
	196	+
	197	+ for (k = 0 ; k < PopSize ; k++)
	198	+ { do
	199	+ a = u_int (double (PopSize) * Rand.UDDouble ()) ;
	200	+ while (a == k) ;
	201	+ do
	202	+ b = u_int (double (PopSize) * Rand.UDDouble ()) ;
	203	+ while (b == k \|\| b == a) ;
	204	+ do
	205	+ c = u_int (double (PopSize) * Rand.UDDouble ()) ;
	206	+ while (c == k \|\| c == a \|\| c == b) ;
	207	+
	208	+if (a >= PopSize) { printf ("a - very bad.\n") ; exit (0) ; } ;
	209	+if (b >= PopSize) { printf ("b - very bad.\n") ; exit (0) ; } ;
	210	+if (c >= PopSize) { printf ("c - very bad.\n") ; exit (0) ; } ;
	211	+
	212	+ scalefactor = 0.5 * (ScaleFactor2 [b] + ScaleFactor2 [c]) ;
	213	+
	214	+ temp = Vector2 [a] + scalefactor * (Vector2 [b] - Vector2 [c]) ;
	215	+
	216	+ crossover = 0.5 * (Crossover2 [k] + Crossover2 [a]) ;
	217	+
	218	+ temp.CrossOver (Vector2 [k], crossover) ;
	219	+
	220	+ cost = Function (temp) ;
	221	+ FuncCalcs ++ ;
	222	+
	223	+ if (cost < Cost2 [k])
	224	+ { MainVector [k] = temp ;
	225	+ MainCost [k] = cost ;
	226	+ ScaleFactor [k] = scalefactor ;
	227	+ Crossover [k] = crossover ;
	228	+ improvement = 0.5 * (improvement + Max (improvement, Cost2 [k] - cost)) ;
	229	+ best = Min (best, cost) ;
	230	+ }
	231	+ else
	232	+ { MainVector [k] = Vector2 [k] ;
	233	+ MainCost [k] = Cost2 [k] ;
	234	+ ScaleFactor [k] = ScaleFactor2 [k] ;
	235	+ Crossover [k] = Crossover2 [k] ;
	236	+ } ; // if
	237	+ } ; // for k
	238	+
	239	+ Generations ++ ;
	240	+ } ; // while
	241	+
	242	+ // Now find best vector.
	243	+
	244	+ for (a = 0, k = 1 ; k < PopSize ; k++)
	245	+ if (MainCost [k] < MainCost [a])
	246	+ a = k ;
	247	+ temp = MainVector [a] ;
	248	+
	249	+ return temp ;
	250	+} ; // DEmin::Minimize
	251	+
	252	+void
	253	+DEmin::PrintStats (void)
	254	+{ printf ("Generations : %8d\n", Generations) ;
	255	+ printf ("Function evaluations : %8d\n", FuncCalcs) ;
	256	+ putchar ('\n') ;
	257	+} ; // DEmin::PrintStats
	258	+
	259	+//======================================================================================
	260	+
	261	+
	262	+
	263	+// Do not edit or modify anything in this comment block.
	264	+// The arch-tag line is a file identity tag for the GNU Arch
	265	+// revision control system.
	266	+//
	267	+// arch-tag: e874ad11-7294-4757-9370-fd1e6e1261d0
	268	+

FIR-POLY/fir_interp.cc+109−0

@@ -0,0 +1,109 @@
	1	+/*
	2	+** Copyright (C) 2004 Erik de Castro Lopo <erikd@mega-nerd.com>
	3	+**
	4	+** The copyright above and this notice must be preserved in all
	5	+** copies of this source code. The copyright above does not
	6	+** evidence any actual or intended publication of this source code.
	7	+**
	8	+** This is unpublished proprietary trade secret of Erik de Castro
	9	+** Lopo. This source code may not be copied, disclosed,
	10	+** distributed, demonstrated or licensed except as authorized by
	11	+** Erik de Castro Lopo.
	12	+**
	13	+** No part of this program or publication may be reproduced,
	14	+** transmitted, transcribed, stored in a retrieval system,
	15	+** or translated into any language or computer language in any
	16	+** form or by any means, electronic, mechanical, magnetic,
	17	+** optical, chemical, manual, or otherwise, without the prior
	18	+** written permission of Erik de Castro Lopo.
	19	+*/
	20	+
	21	+
	22	+#include <cstdio>
	23	+#include <cstdlib>
	24	+#include <ctime>
	25	+#include <cstring>
	26	+#include <complex>
	27	+
	28	+#include <GMatrix.hh>
	29	+#include <Minimize.hh>
	30	+
	31	+#include "mag_spectrum.hh"
	32	+
	33	+
	34	+#define ARRAY_LEN(x) (sizeof (x) / sizeof ((x) [0]))
	35	+
	36	+/*
	37	+** Number of half cycles of impulse response to the left and right
	38	+** of the center.
	39	+** left_half_offset <= right_half_cycles
	40	+*/
	41	+
	42	+// static int left_half_cycles = 5 ;
	43	+// static int right_half_cycles = 15 ;
	44	+
	45	+
	46	+typedef struct
	47	+{
	48	+ double sinc_time_fudge ;
	49	+
	50	+ /* Window coefficients (all in range [0, 1]). */
	51	+ double left_ak [5], right_ak [5] ;
	52	+
	53	+} FIR_PARAMS ;
	54	+
	55	+static double fir_error (const GMatrix& gm) ;
	56	+
	57	+int
	58	+main (void)
	59	+{ static double data [32] ;
	60	+ GMatrix gm_start ;
	61	+ double error ;
	62	+ int k, param_count ;
	63	+
	64	+ param_count = sizeof (FIR_PARAMS) / sizeof (double) ;
	65	+
	66	+ srandom (time (NULL)) ;
	67	+ for (k = 0 ; k < param_count ; k++)
	68	+ data [k] = 3.0 * ((1.0 * random ()) / INT_MAX - 0.5) ;
	69	+
	70	+ gm_start = GMatrix (param_count, 1, data) ;
	71	+
	72	+ fir_error (gm_start) ;
	73	+
	74	+ do
	75	+ error = MinDiffEvol (fir_error, gm_start, 1e-15, 2.0, random ()) ;
	76	+ while (error > 10.0) ;
	77	+
	78	+ printf ("error : %f\n", error) ;
	79	+
	80	+ return 0 ;
	81	+} /* main */
	82	+
	83	+/*==============================================================================
	84	+*/
	85	+
	86	+static double
	87	+fir_error (const GMatrix& gm)
	88	+{ FIR_PARAMS params ;
	89	+ double error = 0.0 ;
	90	+ unsigned param_count ;
	91	+
	92	+
	93	+ param_count = sizeof (params) / sizeof (double) ;
	94	+
	95	+ if (gm.GetData (param_count, (double*) &params) != param_count)
	96	+ { printf ("\n\nError : GetData should return 2.\n") ;
	97	+ exit (1) ;
	98	+ } ;
	99	+
	100	+ return error ;
	101	+} /* fir_error */
	102	+
	103	+/*
	104	+** Do not edit or modify anything in this comment block.
	105	+** The following line is a file identity tag for the GNU Arch
	106	+** revision control system.
	107	+**
	108	+** arch-tag: 3ce7ca6f-394e-432c-aeaa-f228afc05afd
	109	+*/

FIR-POLY/mag_spectrum.cc+62−0

@@ -0,0 +1,62 @@
	1	+/*
	2	+** Copyright (C) 2004 Erik de Castro Lopo <erikd@mega-nerd.com>
	3	+**
	4	+** The copyright above and this notice must be preserved in all
	5	+** copies of this source code. The copyright above does not
	6	+** evidence any actual or intended publication of this source code.
	7	+**
	8	+** This is unpublished proprietary trade secret of Erik de Castro
	9	+** Lopo. This source code may not be copied, disclosed,
	10	+** distributed, demonstrated or licensed except as authorized by
	11	+** Erik de Castro Lopo.
	12	+**
	13	+** No part of this program or publication may be reproduced,
	14	+** transmitted, transcribed, stored in a retrieval system,
	15	+** or translated into any language or computer language in any
	16	+** form or by any means, electronic, mechanical, magnetic,
	17	+** optical, chemical, manual, or otherwise, without the prior
	18	+** written permission of Erik de Castro Lopo.
	19	+*/
	20	+
	21	+#include <cstdio>
	22	+#include <cstdlib>
	23	+#include <cstring>
	24	+#include <cmath>
	25	+
	26	+#include <fftw3.h>
	27	+
	28	+void
	29	+mag_spectrum (double input, int len, double magnitude, double scale)
	30	+{ static fftw_plan plan = NULL ;
	31	+
	32	+ double maxval ;
	33	+ int k ;
	34	+
	35	+ if (input == NULL \|\| magnitude == NULL)
	36	+ return ;
	37	+
	38	+ if (plan == NULL)
	39	+ { plan = fftw_plan_r2r_1d (len, input, magnitude, FFTW_R2HC, FFTW_ESTIMATE \| FFTW_PRESERVE_INPUT) ;
	40	+ if (plan == NULL)
	41	+ { printf ("%s : line %d : create plan failed.\n", __FILE__, __LINE__) ;
	42	+ exit (1) ;
	43	+ } ;
	44	+ } ;
	45	+
	46	+ fftw_execute (plan) ;
	47	+
	48	+ /* (k < N/2 rounded up) */
	49	+ maxval = 0.0 ;
	50	+ for (k = 0 ; k < len / 2 ; k++)
	51	+ magnitude [k] = scale * sqrt (magnitude [k] * magnitude [k] + magnitude [len - k - 1] * magnitude [len - k - 1]) ;
	52	+
	53	+ memset (magnitude + len / 2, 0, len / 2 * sizeof (magnitude [0])) ;
	54	+
	55	+ return ;
	56	+} /* mag_spectrum */
	57	+
	58	+// Do not edit or modify anything in this comment block.
	59	+// The following line is a file identity tag for the GNU Arch
	60	+// revision control system.
	61	+//
	62	+// arch-tag: 58aa43ee-b8e1-4178-8861-7621399fa049

FIR-POLY/mag_spectrum.hh+27−0

@@ -0,0 +1,27 @@
	1	+/*
	2	+** Copyright (C) 2004 Erik de Castro Lopo <erikd@mega-nerd.com>
	3	+**
	4	+** The copyright above and this notice must be preserved in all
	5	+** copies of this source code. The copyright above does not
	6	+** evidence any actual or intended publication of this source code.
	7	+**
	8	+** This is unpublished proprietary trade secret of Erik de Castro
	9	+** Lopo. This source code may not be copied, disclosed,
	10	+** distributed, demonstrated or licensed except as authorized by
	11	+** Erik de Castro Lopo.
	12	+**
	13	+** No part of this program or publication may be reproduced,
	14	+** transmitted, transcribed, stored in a retrieval system,
	15	+** or translated into any language or computer language in any
	16	+** form or by any means, electronic, mechanical, magnetic,
	17	+** optical, chemical, manual, or otherwise, without the prior
	18	+** written permission of Erik de Castro Lopo.
	19	+*/
	20	+
	21	+void mag_spectrum (double input, int len, double magnitude, double scale) ;
	22	+
	23	+// Do not edit or modify anything in this comment block.
	24	+// The following line is a file identity tag for the GNU Arch
	25	+// revision control system.
	26	+//
	27	+// arch-tag: cc5db1c8-4c67-49f7-9c2f-6c5021a69970

Sketch out skeleton of program fir_interp.

Keyboard shortcuts

General

Issues and pull requests

Issues, pull requests and Inbox

Repository code