/ * mkshape -- design raised cosine FIR filter A.J. Fisher, University of York <fisher@minster.york.ac.uk> November 1996 */ #include <stdio.h> #include <math.h> #include "mkfilter.h" #include "complex.h" #define SEQLEN 2048 #define opt_l 0x01 #define opt_c 0x02 /* raised cosine */ #define opt_r 0x04 /* root raised cosine */ #define opt_h 0x08 /* Hilbert transformer */ #define opt_i 0x10 /* identity function */ #define opt_w 0x20 /* apply windowing */ #define opt_x 0x40 /* x / sin x compensation */ #define opt_b 0x80 /* truncate coeffs */ static uint options; static int numcoeffs, truncbits; static double alpha, beta; static double xcoeffs[MAXPZ+1]; static complex circle[SEQLEN/2]; static void readcmdline(char**); static double getf(char*); static int geti(char*); static void usage(), initcircle(), computefilter(), compute_rc(), compute_ht(); static void apply_window(), trunc_coeffs(); static void printresults(char**), printcmdline(char**); static void fft(complex*, complex*, int); static void giveup(char*, int = 0); global void main(int argc, char **argv) { readcmdline(argv); if (options & opt_l) { computefilter(); printresults(argv); } else printf("OK!\n"); exit(0); } static void readcmdline(char **argv) { int ap = 0; unless (argv[ap] == NULL) ap++; /* skip program name */ options = numcoeffs = truncbits = 0; while (argv[ap] != NULL && argv[ap][0] == '-') { char *s = argv[ap++]; s++; uint opts = 0; until (*s == '\0') { char ch = *(s++); if (ch == 'c') opts |= opt_c; else if (ch == 'r') opts |= opt_r; else if (ch == 'h') opts |= opt_h; else if (ch == 'i') opts |= opt_i; else if (ch == 'l') opts |= opt_l; else if (ch == 'x') opts |= opt_x; else if (ch == 'w') opts |= opt_w; else if (ch == 'b') opts |= opt_b; else usage(); } if (opts & (opt_c | opt_r)) { alpha = getf(argv[ap++]); beta = getf(argv[ap++]); } if (opts & (opt_c | opt_r | opt_h | opt_i)) { numcoeffs = geti(argv[ap++]); unless (numcoeffs & 1) giveup("sorry, filter length must be odd"); /* (num.zeros) = (num.coeffs) - 1 */ if (numcoeffs < 1 || numcoeffs > MAXPZ+1) giveup("filter length must be in range 1 .. %d", MAXPZ+1); } if (opts & opt_b) { truncbits = geti(argv[ap++]); if (truncbits <= 0) giveup("Num. bits after -b must be .gt. 0"); } options |= opts; } unless (argv[ap] == NULL) usage(); } static double getf(char *s) { if (s == NULL) usage(); return atof(s); } static int geti(char *s) { if (s == NULL) usage(); return atoi(s); } static void usage() { fprintf(stderr, "Mkshape V.%s from <fisher@minster.york.ac.uk>\n", VERSION); fprintf(stderr, "Usage: mkshape -{cr} <alpha> <beta> <len> [-{lwx}] [-b nb]\n"); fprintf(stderr, " mkshape -i <len> [-{lwx}] [-b nb]\n"); fprintf(stderr, " mkshape -h <len> [-{lw}] [-b nb]\n"); exit(1); } static void initcircle() { for (int i = 0; i < SEQLEN/2; i++) { double x = (TWOPI * i) / SEQLEN; circle[i] = expj(x); /* inverse fft */ } } static void computefilter() { switch (options & (opt_c | opt_r | opt_h | opt_i)) { default: giveup("Must have exactly one of -{crhi}"); case opt_c: case opt_r: case opt_i: initcircle(); compute_rc(); apply_window(); trunc_coeffs(); break; case opt_h: if (options & opt_x) giveup("Can't combine -h with -x"); compute_ht(); apply_window(); trunc_coeffs(); break; } } static void compute_rc() /* (Root?) Raised Cosine */ { complex vec[SEQLEN], temp[SEQLEN]; double f1 = (options & opt_i) ? 0.5 : (1.0-beta)*alpha; double f2 = (options & opt_i) ? 0.5 : (1.0+beta)*alpha; double tau = (options & opt_i) ? 1.0 : 0.5/alpha; for (int i = 0; i <= SEQLEN/2; i++) { double f = (double) i / (double) SEQLEN; vec[i] = (f <= f1) ? 1.0 : (f <= f2) ? 0.5 * (1.0 + cos((PI*tau/beta) * (f-f1))) : 0.0; } if (options & opt_r) { for (int i = 0; i <= SEQLEN/2; i++) vec[i].re = sqrt(vec[i].re); } if (options & opt_x) { for (int i = 1; i <= SEQLEN/2; i++) { double x = PI * (double) i / (double) SEQLEN; vec[i].re *= (x / sin(x)); } } for (int i = 0; i <= SEQLEN/2; i++) vec[i].re *= tau; for (int i = 1; i < SEQLEN/2; i++) vec[SEQLEN-i] = vec[i]; fft(vec, temp, SEQLEN); /* inverse fft */ int h = (numcoeffs-1)/2; for (int i = 0; i < numcoeffs; i++) { int j = (SEQLEN-h+i) % SEQLEN; xcoeffs[i] = vec[j].re / (double) SEQLEN; } } static void compute_ht() /* Hilbert Transformer */ { int h = (numcoeffs-1) / 2; xcoeffs[h] = 0.0; for (int i = 1; i <= h; i++) { double x = (i & 1) ? (1.0 / (double) i) : 0.0; xcoeffs[h+i] = -x; xcoeffs[h-i] = +x; } } static void apply_window() /* apply Hamming window to impulse response */ { if (options & opt_w) { int h = (numcoeffs-1) / 2; for (int i = 1; i <= h; i++) { double w = 0.54 - 0.46 * cos(TWOPI * (double) (h+i) / (double) (numcoeffs-1)); xcoeffs[h+i] *= w; xcoeffs[h-i] *= w; } } } static void trunc_coeffs() { if (options & opt_b) { double fac = pow(2.0, (double) (truncbits-1)); int h = (numcoeffs-1) / 2; double max = (options & opt_h) ? xcoeffs[h-1] : xcoeffs[h]; /* max coeff */ double scale = (fac-1.0) / (fac*max); for (int i = 0; i < numcoeffs; i++) { double x = xcoeffs[i] * scale; /* scale coeffs so max is (fac-1.0)/fac */ xcoeffs[i] = fix(x*fac) / fac; /* truncate */ } } } static void printresults(char **argv) { printcmdline(argv); double gain = 0.0; if (options & opt_h) { int p = ((numcoeffs-1)/2 & 1) ? 0 : 1; bool odd = false; for (int i = p; i < numcoeffs && xcoeffs[i] > 0.0; i += 2) { if (odd) gain += xcoeffs[i]; else gain -= xcoeffs[i]; odd ^= true; } if (odd) gain = -gain; gain *= 2.0; } else { for (int i = 0; i < numcoeffs; i++) gain += xcoeffs[i]; } printf("G = %18.10e\n", gain); printf("NZ = %d\n", numcoeffs-1); for (int i = 0; i < numcoeffs; i++) printf("%18.10e\n", xcoeffs[i]); printf("NP = %d\n", numcoeffs-1); for (int i = 0; i < numcoeffs-1; i++) printf("0\n"); printf("-1\n"); } static void printcmdline(char **argv) { int k = 0; until (argv[k] == NULL) { if (k > 0) putchar(' '); fputs(argv[k++], stdout); } putchar('\n'); } static void fft(complex *data, complex *temp, int n) { if (n > 1) { int h = n/2; for (int i = 0; i < h; i++) { int i2 = i*2; temp[i] = data[i2]; /* even */ temp[h+i] = data[i2+1]; /* odd */ } fft(&temp[0], &data[0], h); fft(&temp[h], &data[h], h); int p = 0, t = SEQLEN/n; for (int i = 0; i < h; i++) { complex wkt = circle[p] * temp[h+i]; data[i] = temp[i] + wkt; data[h+i] = temp[i] - wkt; p += t; } } } static void giveup(char *msg, int p1) { fprintf(stderr, "mkshape: "); fprintf(stderr, msg, p1); putc('\n', stderr); exit(1); }
file: /Techref/uk/ac/york/cs/www-users/http/~fisher/software/mkfilter/current/mkshape.C, 6KB, , updated: 2000/4/4 11:47, local time: 2024/5/1 16:37,
13.58.77.98:LOG IN
|
©2024 These pages are served without commercial sponsorship. (No popup ads, etc...).Bandwidth abuse increases hosting cost forcing sponsorship or shutdown. This server aggressively defends against automated copying for any reason including offline viewing, duplication, etc... Please respect this requirement and DO NOT RIP THIS SITE. Questions? <A HREF="http://www.piclist.com/Techref/uk/ac/york/cs/www-users/http/~fisher/software/mkfilter/current/mkshape.C"> uk ac york cs www-users http ~fisher software mkfilter current mkshape</A> |
Did you find what you needed? |
PICList 2024 contributors:
o List host: MIT, Site host massmind.org, Top posters @none found - Page Editors: James Newton, David Cary, and YOU! * Roman Black of Black Robotics donates from sales of Linistep stepper controller kits. * Ashley Roll of Digital Nemesis donates from sales of RCL-1 RS232 to TTL converters. * Monthly Subscribers: Gregg Rew. on-going support is MOST appreciated! * Contributors: Richard Seriani, Sr. |
Welcome to www.piclist.com! |
.