/* f3.c gcc -Wall -mwindows f3.c f31.o fft/sample1/fft4g.o -lm -lwinmm *** Makefile *** CC = gcc DEFINES = CFLAGS = -Wall -mwindows ${DEFINES} PRGOUT = f3 OBJECTS = f3.o \ f31.o $(PRGOUT): ${OBJECTS} ${CC} $(CFLAGS) -o $(PRGOUT) ${OBJECTS} fft/sample1/fft4g.o -lm -lwinmm .c.o: $(CC) $(CFLAGS) -c $< clean: rm -f *.o $(PRGOUT) */ #include #include #include double itp_val(double*); void send_hLb(char*); void send_logf(char*); int printf_hLb(const char*, ...); /* ---- fft?g functions -------- */ void cdft(int, int, double *, int *, double *); void rdft(int, int, double *, int *, double *); void ddct(int, int, double *, int *, double *); void ddst(int, int, double *, int *, double *); void dfct(int, double *, double *, int *, double *); void dfst(int, double *, double *, int *, double *); /* ---- aio_ functions in f31.c -------- */ int aio_create(HWND); void aio_destroy(); int aio_istart(HWND); int aio_ostart(HWND); int aio_istop(); int aio_ostop(); int aio_WIM_OPEN(); int aio_WIM_DATA(HWND, LPARAM); int aio_WIM_CLOSE(); int aio_WOM_OPEN(); int aio_WOM_DONE(HWND, LPARAM); int aio_WOM_CLOSE(); /* ---- FFT related parameters ------------------ */ /* original value #define NMAX 8192 #define NMAXSQRT 64 */ #define MAX_NPT 16384 #define NMAX (2 * MAX_NPT) #define NMAXSQRT 256 /* 128? */ double a[NMAX + 1]; double b[NMAX + 1]; double w[NMAX * 5 / 4]; int ip[NMAXSQRT + 2] = { 0, }; /* ---- graphic display parameters -------------- */ #define XLEFT 0 /* ƒeƒXƒg—p raw x-scale */ int gp_valid = 0; int gp_npt = 4096; /* gp_npt <= MAX_NPT */ RECT gp; static int gp_width; /* gp.right - gp.left */ static int gp_mh; /* -(gp.bottom - gp.top) */ static int gp_mh2; /* -(gp.bottom - gp.top) / 2 */ static int gp_th2; /* (gp.bottom - gp.top) / 2 + gp.top */ /* ---- make interpolation table --------------- */ #define N_SIDE 3 /* look 3 + 1 + 3 data */ static double itptab[100]; static int itpdiv; /* itp division / 2 */ /* return scale * |z|^2 */ inline double scp(double *p, double scale) { return scale * (*p * *p + *(p + 1) * *(p + 1)); } /* returns scale * |z| */ inline double scabs(double *p, double scale) { return scale * sqrt(*p * *p + *(p + 1) * *(p + 1)); } /* get interpolation value */ inline double itp(double bc) { int i; for (i = 0; i < 2 * itpdiv; i += 1) if (bc < itptab[i]) break; return (double) (i - itpdiv) / (2 * itpdiv); /* 0.01 * (i - 50); */ } inline void add_itp(int f1, double v) { itptab[f1 + itpdiv] += (0.2 * v); /* for 5 frequencies */ } inline double get_f(double f0, int f1) { /* frequency */ return f0 + (0.5 * f1 + 0.25) / itpdiv; } /* make interpolation table */ static int intp(unsigned int d) { if (d > sizeof(itptab)/sizeof(itptab[0])) return -1; itpdiv = d / 2; int f1; /* -itpdiv .. itpdiv - 1 */ double f0; /* normalized fequency */ double v = 4 * 256 * M_PI / gp_npt; /* squere wave amplitude */ int i; double tmp; memset(itptab, 0, sizeof(itptab)); for (f0 = 125; f0 < gp_npt / 2; f0 *= 2) { for (f1 = -itpdiv; f1 < itpdiv; f1 += 1) { #if 1 /* symmetrical signal */ double sft = fmod( get_f(f0, f1) * M_TWOPI * (gp_npt - 1) / gp_npt, M_TWOPI); sft = -sft/2; #else double sft = 0; #endif double* p = a; for (i = 0; i < gp_npt; i += 1) { *p = (cos( get_f(f0, f1) * M_TWOPI * i / gp_npt + sft ) >= 0.0) ? v : -v; p += 2; } /* FFT */ memcpy(b, a, 2 * gp_npt * sizeof(double)); cdft(2 * gp_npt, 1, b, ip, w); /* make table */ double scl = 2.0 / gp_npt; double vmax = 0; int vmax_i = -1; for (i = N_SIDE + 1; i < gp_npt / 2 - N_SIDE; i += 1) { tmp = scp(&b[2 * i], scl); if (tmp > vmax) { vmax = tmp; vmax_i = i; } } double lower = 0, higher = 0; if (vmax_i <= 0) return -1; for (i = vmax_i - N_SIDE; i <= vmax_i + N_SIDE; i += 1) { tmp = scp(&b[2 * i], scl); if (i < vmax_i) lower += tmp; else if (i > vmax_i) higher += tmp; else vmax = tmp; } tmp = (higher - lower) / vmax; add_itp(f1, tmp); // double fin = get_f(f0, f1); // sprintf(cb, "%10.3f, %9g, (%9g, %9g )", // fin, tmp, lower, higher); // send_logf(cb); } // send_logf(""); } #if 0 /* write itp table to log file */ tmp = itptab[0]; for (i = -itpdiv; i < itpdiv; i += 1) { double itptab_i = itptab[i + itpdiv]; sprintf(cb, "%3d: %g (%g, %g) dif %g", i, itptab_i, get_f(0, i), itp(itptab_i), itptab_i - tmp); send_logf(cb); tmp = itptab_i; } send_logf(""); #endif printf_hLb("Interpolation table %d", 2 * itpdiv); return 2 * itpdiv; } /* ---- interpolation ------------------------- */ /* return interpolated value */ double itp_val(double* bp) { double tmp, vmax = 0; int vmax_i = -1; int i; for (i = N_SIDE + 1; i < gp_npt / 2 - N_SIDE; i += 1) { tmp = scp(bp + (2 * i), 1.0); if (tmp > vmax) { vmax = tmp; vmax_i = i; } } if (vmax < 150000 || vmax_i <= 0) return 0; double lower = 0, higher = 0; for (i = vmax_i - N_SIDE; i <= vmax_i + N_SIDE; i += 1) { tmp = scp(bp + (2 * i), 1.0); if (i < vmax_i) lower += tmp; else if (i > vmax_i) higher += tmp; else vmax = tmp; } tmp = (higher - lower) / vmax; return vmax_i + itp(tmp); } /* ---- graphical display ------------------------------ */ /* conversion from double/int to client window cordinate */ inline int cv_ya(double *p, double scale) { return (int)lround(scale * gp_mh * sqrt(*p * *p + *(p + 1) * *(p + 1))) + gp.bottom; } inline int cv_y(double a, double scale) { return (int)lround(scale * gp_mh2 * a) + gp_th2; } inline int cv_x(int i, int n) { #if XLEFT return i + gp.left; #else return (int)(gp_width * i / n) + gp.left; #endif } void on_paint3(HDC hdc) { int i; double* p; POINT pt[gp_npt + 1]; HPEN o_pen; HPEN n_pen1 = CreatePen(PS_DOT, 1, RGB(224, 224, 224)); HPEN n_pen2 = CreatePen(PS_SOLID, 1, RGB(224, 224, 224)); HPEN n_pen3 = CreatePen(PS_SOLID, 1, RGB(0, 0, 255)); #if 1 /* imaginary part */ p = a + 1; for (i = 0; i < gp_npt; i += 1) { pt[i].x = cv_x(i, gp_npt); pt[i].y = cv_y(*p, 1.0); p += 2; } pt[i].x = pt[i - 1].x + 1; pt[i].y = pt[i - 1].y; o_pen = (HPEN) SelectObject(hdc, n_pen1); Polyline(hdc, pt, gp_npt + 1); /* real part */ p = a; for (i = 0; i < gp_npt; i += 1) { /* pt[i].x = cv_x(i, gp_npt); */ pt[i].y = cv_y(*p, 1.0); p += 2; } /* pt[i].x = pt[i - 1].x + 1; */ pt[i].y = pt[i - 1].y; SelectObject(hdc, n_pen2); Polyline(hdc, pt, gp_npt + 1); #endif /* amplitude spectrum */ p = b; SelectObject(hdc, n_pen3); for (i = 0; i < gp_npt / 2; i += 1) { pt[0].x = cv_x(i, gp_npt / 2); pt[0].y = gp.bottom; pt[1].x = pt[0].x; pt[1].y = cv_ya(p, 2.0 / gp_npt); pt[2].x = pt[1].x; pt[2].y = pt[1].y + 1; Polyline(hdc, pt, 3); p += 2; } SelectObject(hdc, o_pen); DeleteObject(n_pen3); DeleteObject(n_pen2); DeleteObject(n_pen1); } void on_paint2(HDC hdc) { int i; POINT pt[4], *p; HPEN n_pen1 = CreatePen(PS_DOT, 1, RGB(192, 192, 192)); HPEN n_pen2 = CreatePen(PS_SOLID, 1, RGB(192, 192, 192)); HPEN o_pen = (HPEN) SelectObject(hdc, n_pen1); RECT rt = gp; /* vertical center line */ p = (POINT *) & pt; i = (rt.left + rt.right) / 2; p->x = i, p->y = rt.bottom, p += 1; p->x = i, p->y = rt.top, p += 1; Polyline(hdc, pt, p - pt); /* horizontal center line */ p = (POINT *) & pt; i = (rt.top + rt.bottom) / 2; p->x = rt.left, p->y = i, p += 1; p->x = rt.right, p->y = i, p += 1; Polyline(hdc, pt, p - pt); /* top and right line (dotted) */ p = (POINT *) & pt; p->x = rt.left, p->y = rt.top, p += 1; p->x = rt.right, p->y = rt.top, p += 1; p->x = rt.right, p->y = rt.bottom, p += 1; Polyline(hdc, pt, p - pt); /* left and bottom line */ SelectObject(hdc, n_pen2); p = (POINT *) & pt; p->x = rt.left, p->y = rt.top, p += 1; p->x = rt.left, p->y = rt.bottom, p += 1; p->x = rt.right, p->y = rt.bottom, p += 1; p->x = (p - 1)->x + 1, p->y = (p - 1)->y, p += 1; /* added */ Polyline(hdc, pt, p - pt); /* graduation of scale */ SelectObject(hdc, o_pen); for (i = rt.bottom; i >= rt.top; i += gp_mh / 10) { /* SetPixel(hdc, rt.left - 1, i, RGB(0, 0, 0)); SetPixel(hdc, rt.left - 2, i, RGB(0, 0, 0)); */ p = (POINT *) & pt; p->x = rt.left - 1, p->y = i, p += 1; p->x = rt.left - 3, p->y = i, p += 1; Polyline(hdc, pt, p - pt); } for (i = rt.left; i <= rt.right; i += gp_width / 16) { /* SetPixel(hdc, i, rt.bottom + 1, RGB(0, 0, 0)); SetPixel(hdc, i, rt.bottom + 2, RGB(0, 0, 0)); */ p = (POINT *) & pt; p->x = i, p->y = rt.bottom + 1, p += 1; p->x = i, p->y = rt.bottom + 3, p += 1; Polyline(hdc, pt, p - pt); } DeleteObject(n_pen2); DeleteObject(n_pen1); } void on_paint(HDC hdc, RECT* rp) { if (gp_valid) { on_paint2(hdc); on_paint3(hdc); } } /* ---- other window message processing --------------- */ #define TID(n) (n) static HFONT fnt; static HWND hLb; void on_timer(HWND hWnd, int id) { if (id == 1) { /* from WM_CREATE */ send_hLb("FFTA"); RECT rt; GetClientRect(hWnd, &rt); gp.left = rt.left + 60; gp.bottom = rt.bottom - 40; gp.right = gp.left + 512; gp.top = gp.bottom - 400; gp_width = gp.right - gp.left; gp_mh = -(gp.bottom - gp.top); gp_mh2 = gp_mh / 2; gp_th2 = (gp.bottom - gp.top) / 2 + gp.top; gp_valid = 1; SetCursor(LoadCursor(NULL, IDC_WAIT)); SetTimer(hWnd, TID(2), 50, NULL); InvalidateRect(hWnd, NULL, TRUE); } else if (id == 2) { intp(100); memset(a, 0, (2 * gp_npt + 1) * sizeof(double)); memset(b, 0, (2 * gp_npt + 1) * sizeof(double)); SetCursor(LoadCursor(NULL, IDC_ARROW)); SetTimer(hWnd, TID(3), 100, NULL); InvalidateRect(hWnd, NULL, FALSE); } else if (id == 3) { aio_istart(hWnd); } } int on_create(HWND hWnd, WPARAM wParam, LPARAM lParam) { LOGFONT lf; lf.lfHeight = 16; lf.lfWidth = 0; lf.lfEscapement = lf.lfOrientation = 0; lf.lfWeight = FW_NORMAL; lf.lfItalic = lf.lfUnderline = FALSE; lf.lfStrikeOut = FALSE; lf.lfCharSet = DEFAULT_CHARSET; lf.lfOutPrecision = OUT_DEFAULT_PRECIS; lf.lfClipPrecision = CLIP_DEFAULT_PRECIS; lf.lfQuality = DEFAULT_QUALITY; lf.lfPitchAndFamily = FIXED_PITCH; lf.lfFaceName[0] = '\0'; fnt = CreateFontIndirect(&lf); hLb = CreateWindow(TEXT("LISTBOX"), NULL, WS_CHILD /* | WS_POPUP */ | WS_VISIBLE | WS_VSCROLL | WS_HSCROLL | WS_BORDER | WS_CAPTION | LBS_NOTIFY, 900 - 320, 1, 480, 320, hWnd, (HMENU) 1, ((LPCREATESTRUCT) (lParam))->hInstance, NULL); SendMessage(hLb, WM_SETFONT, (WPARAM) fnt, FALSE); SendMessage(hLb, LB_SETHORIZONTALEXTENT, 3 * 320, 0); SetTimer(hWnd, TID(1), 50, NULL); if (!hLb) return -1; return aio_create(hWnd); } void on_char(HWND hWnd, int id) { } /* ---- report functions --------------------------------- */ void send_hLb(char* s) { SendMessage(hLb, LB_ADDSTRING, 0, (LPARAM) s); PostMessage(hLb, LB_SETTOPINDEX, (WPARAM) (SendMessage(hLb, LB_GETCOUNT, 0, 0) - 1), 0); } void send_logf(char* s) { static FILE* of = NULL; if (!of) of = fopen("log", "w"); if (!of) return; fputs(s, of); fputc('\n', of); } int printf_hLb(const char* fmt, ...) { char cb[4096]; va_list ap; va_start(ap, fmt); int r = vsnprintf(cb, sizeof(cb), fmt, ap); va_end(ap); send_hLb(cb); return r; } /* ---- WndProc and WinMain ----------------------------- */ static char glpszText[256]; LRESULT CALLBACK WndProc(HWND hWnd, UINT message, WPARAM wParam, LPARAM lParam) { PAINTSTRUCT ps; HDC hdc; RECT rt; switch (message) { case MM_WIM_OPEN: aio_WIM_OPEN(hWnd); return 0; case MM_WIM_DATA: aio_WIM_DATA(hWnd, lParam); return 0; case MM_WIM_CLOSE: aio_WIM_CLOSE(); return 0; case MM_WOM_OPEN: aio_WOM_OPEN(hWnd); return 0; case MM_WOM_DONE: aio_WOM_DONE(hWnd, lParam); return 0; case MM_WOM_CLOSE: aio_WOM_CLOSE(); return 0; case WM_TIMER: KillTimer(hWnd, wParam); if (TID(1) <= wParam && wParam <= TID(3)) { on_timer(hWnd, wParam); return 0; } break; case WM_PAINT: hdc = BeginPaint(hWnd, &ps); GetClientRect(hWnd, &rt); DrawText(hdc, glpszText, strlen(glpszText), &rt, DT_TOP | DT_LEFT); on_paint(hdc, &rt); EndPaint(hWnd, &ps); return 0; case WM_CHAR: on_char(hWnd, wParam); break; case WM_KEYDOWN: if (wParam == VK_ESCAPE) { PostMessage(hWnd, WM_CLOSE, 0, 0); return 0; } break; case WM_DESTROY: aio_istop(); aio_destroy(); DeleteObject(fnt); PostQuitMessage(0); return 0; case WM_CREATE: on_create(hWnd, wParam, lParam); /* break; */ } return DefWindowProc(hWnd, message, wParam, lParam); } int APIENTRY WinMain(HINSTANCE hInstance, HINSTANCE hPrevInstance, LPSTR lpCmdLine, int nCmdShow) { sprintf(glpszText, "GetCommandLine(): [%s]\n" "WinMain lpCmdLine: [%s]\n" "MAX_NPT %d. Press ESC to exit.\n", lpCmdLine, GetCommandLine(), MAX_NPT); WNDCLASSEX wcex; wcex.cbSize = sizeof(wcex); wcex.style = CS_HREDRAW | CS_VREDRAW; wcex.lpfnWndProc = WndProc; wcex.cbClsExtra = 0; wcex.cbWndExtra = 0; wcex.hInstance = hInstance; wcex.hIcon = LoadIcon(NULL, IDI_APPLICATION); wcex.hCursor = LoadCursor(NULL, IDC_ARROW); wcex.hbrBackground = (HBRUSH) (COLOR_WINDOW + 1); wcex.lpszMenuName = NULL; wcex.lpszClassName = "FFTA"; wcex.hIconSm = NULL; if (!RegisterClassEx(&wcex)) return FALSE; HWND hWnd = CreateWindow("FFTA", "FFTA", WS_OVERLAPPEDWINDOW, CW_USEDEFAULT, CW_USEDEFAULT, CW_USEDEFAULT, CW_USEDEFAULT, NULL, NULL, hInstance, NULL); if (!hWnd) return FALSE; ShowWindow(hWnd, nCmdShow); UpdateWindow(hWnd); MSG msg; while (GetMessage(&msg, NULL, 0, 0)) { TranslateMessage(&msg); DispatchMessage(&msg); } return msg.wParam; }