home *** CD-ROM | disk | FTP | other *** search
/ The Datafile PD-CD 4 / DATAFILE_PDCD4.iso / utilities / utilsf / fftgraph / !FFT_2D / scr / c / FFT_2D
Encoding:
Text File  |  1995-09-16  |  14.9 KB  |  508 lines
  1. /* fft_frequ.c */
  2.  
  3. #include <stdio.h>
  4. #include <stdlib.h>
  5. #include <math.h>
  6. #include <string.h>
  7. #include <kernel.h>
  8. #include <swis.h>
  9. /* #include "visdelay.h" */
  10. #include "fft.h"
  11.  
  12.  
  13. typedef struct {
  14.   char name[12];
  15.   float rate;
  16.   unsigned int window;
  17.   unsigned int nr_of_windows;
  18. } file_header;
  19.  
  20. typedef struct {
  21.   unsigned int draw;
  22.   unsigned int version;
  23.   unsigned int subversion;
  24.   char creator[12];
  25.   int low_x;
  26.   int low_y;
  27.   int hi_x;
  28.   int hi_y;
  29.   int misc[22];
  30. } DrawHeader;
  31.  
  32. extern MCPower(float *a, float *b, int count);
  33. extern MCHanning(float *hanning, int n_data);
  34. int MakDraw(file_header fheader, float *in_data, char *draw_name, float s);
  35.  
  36. int main(int argc, char *argv[]);
  37.  
  38. int main(int argc, char *argv[])
  39. {
  40.  
  41.   float *real, *imag, *hanning, rate, scaleoffset;
  42.   int m, n_data, tmp, type, offset;
  43.   register int n;
  44.   long int flen, ltmp;
  45.   signed char schar;
  46.   FILE *in;
  47.   file_header header;
  48.   _kernel_swi_regs reg;
  49.   _kernel_oserror *err;
  50.  
  51.  
  52.   if(argc != 5) {
  53.     fprintf(stderr, "Usage: SampleFile DrawFile ZeroOffset ScalingFactor\n");
  54.     exit(EXIT_FAILURE);
  55.   }
  56.   if((in = fopen(argv[1], "rb")) == NULL) {
  57.     fprintf(stderr, "Cannot open %s for input\n", argv[1]);
  58.     exit(EXIT_FAILURE);
  59.   }
  60.  
  61.   reg.r[0] = 17;
  62.   reg.r[1] = (int)argv[1];
  63.   if((err=_kernel_swi(OS_File,®,®))!=NULL) {
  64.     fprintf(stderr,"%s\n",err->errmess);
  65.     exit(EXIT_FAILURE);
  66.   }
  67.  
  68.   flen = reg.r[4];
  69.   type = (reg.r[2] & 0xFFF00) >> 8;
  70.  
  71.   if(type != 0xD3C) {
  72.     fprintf(stderr,"Input file must be an Armadeus one (0xD3C)\n");
  73.     exit(EXIT_FAILURE);
  74.   }
  75.   offset = atoi(argv[3]);
  76.   scaleoffset = atof(argv[4]);
  77.   if(scaleoffset == 0) {
  78.     fprintf(stderr,"You can't scale by 0, do you ?\n");
  79.     exit(EXIT_FAILURE);
  80.   }
  81.  
  82.   printf("Input : %s\nType  : %X (Armadeus)\nLen   : %ld\n\n",argv[1],type,flen);
  83.   printf("Output: %s\nType  : AFF (Draw)\n\n",argv[2]);
  84.  
  85.   tmp = 1;
  86.   ltmp = flen-2;
  87.   while(ltmp >>= 1) tmp++;
  88.   m = tmp;
  89.   n_data = 1;
  90.   while(tmp--) n_data <<= 1;
  91.  
  92. /*  visdelay_init(); */
  93.  
  94.   if((real = malloc(sizeof(float) * n_data)) == NULL) {
  95.     fprintf(stderr,"Cannot allocate %d bytes for real\n", sizeof(float) * n_data);
  96.     exit(EXIT_FAILURE);
  97.   }
  98.   if((imag = malloc(sizeof(float) * n_data)) == NULL) {
  99.     fprintf(stderr,"Cannot allocate %d bytes for imag\n", sizeof(float) * n_data);
  100.     exit(EXIT_FAILURE);
  101.   }
  102.   if((hanning = malloc(sizeof(float) * n_data)) == NULL) {
  103.     fprintf(stderr,"Cannot allocate %d bytes for hanning\n", sizeof(float) * n_data);
  104.     exit(EXIT_FAILURE);
  105.   }
  106.  
  107. /*  visdelay_begin(); */
  108.  
  109.   for(n=0; n<n_data; n++) {
  110.     real[n] = 0;
  111.     imag[n] = 0;
  112.   }
  113.  
  114.   MCHanning(hanning, n_data);
  115.  
  116.   header.rate = (float)getc(in) * (float)1E-6;
  117.   printf("Windowing                ");
  118.   if(offset != 0) {
  119.      for(n=0; n<flen-1; n++) {
  120.         schar = (signed char)getc(in);
  121.         real[n] = (float)((int)schar + offset) * hanning[n];
  122.      }
  123.   }
  124.   else {
  125.      for(n=0; n<flen-1; n++) {
  126.         schar = (signed char)getc(in);
  127.         real[n] = (float)(schar) * hanning[n];
  128.      }
  129.   }
  130.   fclose(in);
  131.   printf("- Done\nComputing FFT            ");
  132.   fft2(real, imag, m-1);
  133.   printf("- Done\nComputing power spectrum ");
  134.   MCPower(real, imag, n_data>>1);
  135.   printf("- Done\n");
  136.   strcpy(header.name, "FFT_Frequ");
  137. /*  header.rate = rate; */
  138.   header.window = n_data >> 1;
  139.   header.nr_of_windows = 1;
  140.  
  141.   free(hanning);
  142.   free(imag);
  143. /*  fwrite(&header, sizeof(file_header), 1, out); */
  144. /*  fwrite(real, sizeof(float), n_data>>1, tmp);  */
  145.   printf("Building graph           ");
  146.   MakDraw(header, real, argv[2], scaleoffset);
  147.   printf("- Done\n");
  148. /*  fclose(out); */
  149.  
  150. /*  visdelay_end(); */
  151.  
  152.   return(0);
  153.   }
  154.  
  155. /*****************************************************************************/
  156.  
  157. int MakDraw(file_header fheader, float *in_data, char *draw_name, float s)
  158. {
  159.   FILE *out;
  160.   DrawHeader header;
  161.   unsigned int *buffer, buff_len;
  162.   unsigned int paper_height, paper_width, usable_height, usable_width;
  163.   unsigned int bottom_margin, left_margin, height, mm, point;
  164.   unsigned int flen, ptr, start, end, tag_nr, tag_len;
  165.   unsigned int x, x_origin, y, y_origin;
  166.   int n, m;
  167.  
  168.   char tag_label[80], date[40];
  169.   float  y_max;
  170.   double low_x, high_x, x_factor, y_factor, tag_step, nyquist;
  171.  
  172.   _kernel_swi_regs reg;
  173.   _kernel_oserror *err;
  174.  
  175.  
  176.  
  177.   if((out=fopen(draw_name,"wb"))==NULL) {
  178.     fprintf(stderr,"Cannot open %s for output\n", draw_name);
  179.     exit(EXIT_FAILURE);
  180.   }
  181.  
  182.   flen = fheader.window * fheader.nr_of_windows;
  183.  
  184.   buff_len = 3 * fheader.window * sizeof(int) + 51200;
  185.   if((buffer=malloc(buff_len))==NULL) {
  186.     fprintf(stderr,"Cannot allocate %d bytes for buffer\n",buff_len);
  187.     exit(EXIT_FAILURE);
  188.   }
  189.   nyquist = 0.5 / (double)fheader.rate;
  190.  
  191.   header.draw       = 0x77617244;
  192.   header.version    = 0x000000C9;
  193.   header.subversion = 0x00000000;
  194.   strcpy(header.creator,fheader.name);
  195.   header.low_x      = 0x00000000;
  196.   header.low_y      = 0x00000000;
  197.   header.hi_x       = 8 * 0xB400;
  198.   header.hi_y       = (unsigned int)(11.6 * 0xB400);
  199.   header.misc[0]    = 0x0000000B;        /* object type 11 */
  200.   header.misc[1]    = 0x00000058;    /* object size */
  201.   header.misc[2]    = 0x00000000;
  202.   header.misc[3]    = 0x00000000;
  203.   header.misc[4]    = 0x00000000;
  204.   header.misc[5]    = 0x00000000;
  205.   header.misc[6]    = 0x00000500;    /* paper size */
  206.   header.misc[7]    = 0x00000101;    /* grid default + portrait + limit shown */
  207.   header.misc[8]    = 0x3FF00000;    /* gid spacing */
  208.   header.misc[9]    = 0x00000000;
  209.   header.misc[10]   = 0x00000002;    /* grid division */
  210.   header.misc[11]   = 0x00000000;    /* grid rectangular */
  211.   header.misc[12]   = 0x00000000;    /* grid no auto adjust */
  212.   header.misc[13]   = 0x00000000;    /* grid no show */
  213.   header.misc[14]   = 0x00000000;    /* grid no lock */
  214.   header.misc[15]   = 0x00000001;    /* grid cm */
  215.   header.misc[16]   = 0x00000001;    /* zoom ratio numerator */
  216.   header.misc[17]   = 0x00000001;    /* zoom ratio denominator */
  217.   header.misc[18]   = 0x00000000;    /* no lock */
  218.   header.misc[19]   = 0x00000001;    /* toolbox present */
  219.   header.misc[20]   = 0x00000080;    /* select pre-selected */
  220.   header.misc[21]   = 0x00001388;    /* undo buffer size */
  221.   fwrite(&header, sizeof(DrawHeader), 1, out);
  222.  
  223.   ptr = 0;
  224.   paper_width = 0x5A000;
  225.   paper_height= 0x82800;
  226.   height = (paper_width * 3) >> 2;
  227.   bottom_margin = 0x5A00;
  228.   left_margin = 0x16AE;
  229.   mm = 1814;
  230.   point = 640;
  231.   y_max = 0;
  232.   x_origin = left_margin + 10*mm;
  233.   y_origin = bottom_margin + 10*mm;
  234.   usable_width = paper_width - 20*mm - (left_margin<<1);
  235.   usable_height = height - (y_origin<<1);
  236.  
  237. #ifdef DEBUG
  238.    fprintf(stderr,"usable_width: %u\n",usable_width);
  239. #endif
  240.  
  241. /* horizontal line */
  242.   start = ptr;
  243.   buffer[ptr++] = 0x00000002;                        /* Path object            */
  244.   buffer[ptr++] = 0x00000044;                        /* object size            */
  245.   buffer[ptr++] = x_origin;                        /* bounding box low_x        */
  246.   buffer[ptr++] = y_origin;                        /*              low_y        */
  247.   buffer[ptr++] = x_origin + usable_width;                /*        high_x        */
  248.   buffer[ptr++] = y_origin + 0;                        /*        high_y        */
  249.   buffer[ptr++] = 0xFFFFFFFF;                        /* no fill            */
  250.   buffer[ptr++] = 0x00000000;                        /* outline colour        */
  251.   buffer[ptr++] = point >> 2;                        /* outline width 1/4 point    */
  252.   buffer[ptr++] = 0x20100042;                        /* style            */
  253.   buffer[ptr++] = 0x00000002;                        /* move to            */
  254.   buffer[ptr++] = x_origin;                        /*       x            */
  255.   buffer[ptr++] = y_origin;                        /*       y            */
  256.   buffer[ptr++] = 0x00000008;                        /* line to            */
  257.   buffer[ptr++] = x_origin + usable_width;                /*       x            */
  258.   buffer[ptr++] = y_origin + 0;                        /*       y            */
  259.   buffer[ptr++] = 0x00000000;                        /* end path            */
  260.   end = ptr;
  261.   buffer[start+1] = (end-start) << 2;
  262.  
  263.   fwrite(buffer, sizeof(int), ptr, out);
  264.   ptr=0;
  265.  
  266. /* vertical line */
  267.   start = ptr;
  268.   buffer[ptr++] = 0x00000002;                        /* Path object            */
  269.   buffer[ptr++] = 0x00000000;                        /* object size            */
  270.   buffer[ptr++] = x_origin;                        /* bounding box low_x        */
  271.   buffer[ptr++] = y_origin;                        /*              low_y        */
  272.   buffer[ptr++] = x_origin + 0;                        /*        high_x        */
  273.   buffer[ptr++] = y_origin + usable_height;                /*        high_y        */
  274.   buffer[ptr++] = 0xFFFFFFFF;                        /* no fill            */
  275.   buffer[ptr++] = 0x00000000;                        /* outline colour        */
  276.   buffer[ptr++] = point >> 2;                        /* outline width 1/4 point    */
  277.   buffer[ptr++] = 0x20100042;                        /* style            */
  278.   buffer[ptr++] = 0x00000002;                        /* move to            */
  279.   buffer[ptr++] = x_origin;                        /*       x            */
  280.   buffer[ptr++] = y_origin;                        /*       y            */
  281.   buffer[ptr++] = 0x00000008;                        /* line to            */
  282.   buffer[ptr++] = x_origin + 0;                        /*       x            */
  283.   buffer[ptr++] = y_origin + usable_height;                /*       y            */
  284.   buffer[ptr++] = 0x00000000;                        /* end path            */
  285.   end = ptr;
  286.   buffer[start+1] = (end-start) << 2;
  287.  
  288.   fwrite(buffer, sizeof(int), ptr, out);
  289.   ptr=0;
  290.  
  291. /* filename */
  292.   reg.r[0] = 14;
  293.   reg.r[1] = (int)date;
  294.   date[0] = 0;
  295.   if((err=_kernel_swi(OS_Word, ®, ®)) != NULL) {
  296.     fprintf(stderr,"%s\n", err->errmess);
  297.     exit(EXIT_FAILURE);}
  298.  
  299.   date[24] = ')';
  300.   date[25] = '\0';
  301.  
  302.   sprintf(tag_label,"Name : %s    (",draw_name);
  303.   strcat(tag_label,date);
  304.   tag_len = strlen(tag_label)+1;
  305.   start = ptr;
  306.   buffer[ptr++] = 0x00000001;
  307.   buffer[ptr++] = 0x00000000;
  308.   buffer[ptr++] = x_origin;
  309.   buffer[ptr++] = height - 8*mm - 0x2000;
  310.   buffer[ptr++] = x_origin + tag_len*0x1000;
  311.   buffer[ptr++] = height - 8*mm;
  312.   buffer[ptr++] = 0x00000000;
  313.   buffer[ptr++] = 0xFFFFFF00;
  314.   buffer[ptr++] = 0x00000000;
  315.   buffer[ptr++] = 0x00001000;
  316.   buffer[ptr++] = 0x00002000;
  317.   buffer[ptr++] = x_origin;
  318.   buffer[ptr++] = height - 8*mm - 0x1C00;
  319.   for(m=0; m<=tag_len/4; m++) buffer[ptr+m] = 0;
  320.   strcpy((char *)(buffer+ptr),tag_label);
  321.  
  322.   ptr += tag_len / 4 + 1;
  323.   end = ptr;
  324.   buffer[start+1] = (end-start) << 2;
  325.  
  326.   fwrite(buffer, sizeof(int), ptr, out);
  327.   ptr=0;
  328.  
  329.  
  330. /* 100 Hz tags */
  331.   tag_step = (100.0 * (double)usable_width) / nyquist;
  332.   tag_nr = (unsigned int)((double)usable_width / tag_step);
  333.  
  334.   for(n=1; n<=tag_nr; n++) {
  335.  
  336.     low_x = tag_step * n;
  337.     high_x = low_x;
  338.  
  339.     start = ptr;
  340.     buffer[ptr++] = 0x00000002;
  341.     buffer[ptr++] = 0x00000000;
  342.     buffer[ptr++] = x_origin + (unsigned int)low_x;
  343.     buffer[ptr++] = y_origin - 1*mm;
  344.     buffer[ptr++] = x_origin + (unsigned int)high_x;
  345.     buffer[ptr++] = y_origin;
  346.     buffer[ptr++] = 0xFFFFFFFF;
  347.     buffer[ptr++] = 0x00000000;
  348.     buffer[ptr++] = point >> 2;
  349.     buffer[ptr++] = 0x20100042;
  350.     buffer[ptr++] = 0x00000002;
  351.     buffer[ptr++] = x_origin + (unsigned int)low_x;
  352.     buffer[ptr++] = y_origin - 1*mm;
  353.     buffer[ptr++] = 0x00000008;
  354.     buffer[ptr++] = x_origin + (unsigned int)high_x;
  355.     buffer[ptr++] = y_origin;
  356.     buffer[ptr++] = 0x00000000;
  357.     end = ptr;
  358.     buffer[start+1] = (end-start) << 2;
  359.  
  360.     fwrite(buffer, sizeof(int), ptr, out);
  361.     ptr=0;
  362.   }
  363.  
  364. /* 1000 Hz tags */
  365.   tag_step = (1000.0 * (double)usable_width) / nyquist;
  366.   tag_nr = (unsigned int)((double)usable_width / tag_step);
  367.  
  368.   for(n=0; n<=tag_nr; n++) {
  369.     low_x = tag_step * n;
  370.     high_x = low_x;
  371.  
  372.     start = ptr;
  373.     buffer[ptr++] = 0x00000002;
  374.     buffer[ptr++] = 0x00000000;
  375.     buffer[ptr++] = x_origin + (unsigned int)low_x;
  376.     buffer[ptr++] = y_origin - 3*mm;
  377.     buffer[ptr++] = x_origin + (unsigned int)high_x;
  378.     buffer[ptr++] = y_origin;
  379.     buffer[ptr++] = 0xFFFFFFFF;
  380.     buffer[ptr++] = 0x00000000;
  381.     buffer[ptr++] = point >> 2;
  382.     buffer[ptr++] = 0x20100042;
  383.     buffer[ptr++] = 0x00000002;
  384.     buffer[ptr++] = x_origin + (unsigned int)low_x;
  385.     buffer[ptr++] = y_origin - 3*mm;
  386.     buffer[ptr++] = 0x00000008;
  387.     buffer[ptr++] = x_origin + (unsigned int)high_x;
  388.     buffer[ptr++] = y_origin;
  389.     buffer[ptr++] = 0x00000000;
  390.     end = ptr;
  391.     buffer[start+1] = (end-start) << 2;
  392.  
  393.     fwrite(buffer, sizeof(int), ptr, out);
  394.     ptr=0;
  395.   }
  396.  
  397. /* tag_label */
  398.   for(n=0; n<=tag_nr; n++) {
  399.     sprintf(tag_label,"%d",(int)((float)n*s));
  400.     tag_len = strlen(tag_label);
  401.     low_x = tag_step * n;
  402.     high_x = low_x;
  403.  
  404.     start = ptr;
  405.     buffer[ptr++] = 0x00000001;
  406.     buffer[ptr++] = 0x00000000;
  407.     buffer[ptr++] = x_origin + (unsigned int)low_x - tag_len*0x800;
  408.     buffer[ptr++] = y_origin - 4*mm - 0x2000;
  409.     buffer[ptr++] = x_origin + (unsigned int)high_x + tag_len*0x800;
  410.     buffer[ptr++] = y_origin - 4*mm;
  411.     buffer[ptr++] = 0x00000000;
  412.     buffer[ptr++] = 0xFFFFFF00;
  413.     buffer[ptr++] = 0x00000000;
  414.     buffer[ptr++] = 0x00001000;
  415.     buffer[ptr++] = 0x00001000;
  416.     buffer[ptr++] = x_origin + (unsigned int)low_x - tag_len*0x800;
  417.     buffer[ptr++] = y_origin - 4*mm - 0x1C00;
  418.     for(m=0; m<=tag_len>>2; m++) buffer[ptr+m] = 0;
  419.     strcpy((char *)(buffer+ptr),tag_label);
  420.     ptr += (tag_len>>2) + 1;
  421.     end = ptr;
  422.     buffer[start+1] = (end-start) << 2;
  423.  
  424.     fwrite(buffer, sizeof(int), ptr, out);
  425.     ptr=0;
  426.   }
  427.  
  428. /* kHz */
  429.   sprintf(tag_label,"kHz");
  430.   tag_len = strlen(tag_label);
  431.   low_x = tag_step * (tag_nr+1.0);
  432.   high_x = low_x;
  433.  
  434.   start = ptr;
  435.   buffer[ptr++] = 0x00000001;
  436.   buffer[ptr++] = 0x00000000;
  437.   buffer[ptr++] = x_origin + (unsigned int)low_x - tag_len*0x800;
  438.   buffer[ptr++] = y_origin - 4*mm - 0x2000;
  439.   buffer[ptr++] = x_origin + (unsigned int)high_x + tag_len*0x800;
  440.   buffer[ptr++] = y_origin - 4*mm;
  441.   buffer[ptr++] = 0x00000000;
  442.   buffer[ptr++] = 0xFFFFFF00;
  443.   buffer[ptr++] = 0x00000000;
  444.   buffer[ptr++] = 0x00001000;
  445.   buffer[ptr++] = 0x00001000;
  446.   buffer[ptr++] = x_origin + (unsigned int)low_x - tag_len*0x800;
  447.   buffer[ptr++] = y_origin - 4*mm - 0x1C00;
  448.   for(m=0; m<=tag_len>>2; m++) buffer[ptr+m] = 0;
  449.   strcpy((char *)(buffer+ptr),tag_label);
  450.   ptr += (tag_len>>2) + 1;
  451.   end = ptr;
  452.   buffer[start+1] = (end-start) << 2;
  453.  
  454.   fwrite(buffer, sizeof(int), ptr, out);
  455.   ptr=0;
  456.  
  457.   for(n=1; n<flen; n++) if(in_data[n] > y_max) y_max = in_data[n];
  458.   y_factor = (double)usable_height / (double)y_max;
  459.   x_factor = (double)usable_width / (double)flen;
  460.  
  461. /* diagram */
  462.   start = ptr;
  463.   buffer[ptr++] = 0x00000002;                        /* Path object            */
  464.   buffer[ptr++] = 0x00000000;                        /* object size            */
  465.   buffer[ptr++] = x_origin;                        /* bounding box low_x        */
  466.   buffer[ptr++] = y_origin;                        /*              low_y        */
  467.   buffer[ptr++] = x_origin;                        /*        high_x        */
  468.   buffer[ptr++] = y_origin;                        /*        high_y        */
  469.   buffer[ptr++] = 0xFFFFFFFF;                        /* no fill            */
  470.   buffer[ptr++] = 0x00000000;                        /* outline colour        */
  471.   buffer[ptr++] = 0x00000000; /*point >> 2;*/                /* outline width 1/4 point    */
  472.   buffer[ptr++] = 0x20100042;                        /* style            */
  473.   buffer[ptr++] = 0x00000002;                        /* move to            */
  474.   buffer[ptr++] = x_origin;                        /*       x            */
  475.   buffer[ptr++] = y_origin;                        /*       y            */
  476.  
  477.   for(n=1; n<flen;n++) {
  478.     buffer[ptr++] = 0x00000008;                        /* line to            */
  479.     x = x_origin + (unsigned int)((n-1.0)*x_factor);
  480.     if(x>x_origin) buffer[start+4] = x;
  481.     buffer[ptr++] = x;                            /*       x            */
  482.     y = (int)((double)in_data[n] * y_factor) + y_origin;
  483.     if(y>y_origin) buffer[start+5] = y;
  484.     buffer[ptr++] = y;                            /*       y            */
  485.   }
  486.  
  487.   buffer[ptr++] = 0x00000000;
  488.   end = ptr;
  489.   buffer[start+1] = (end-start) << 2;
  490.  
  491.  
  492.  
  493.   fwrite(buffer, sizeof(int), ptr, out);
  494.   fclose(out);
  495.  
  496.   reg.r[0] = 18;
  497.   reg.r[1] = (int)draw_name;
  498.   reg.r[2] = 0xAFF;
  499.   if((err=_kernel_swi(OS_File,®,®))!=NULL) {
  500.     fprintf(stderr,"Error : %s\n",err->errmess);
  501.     exit(EXIT_FAILURE);
  502.   }
  503. /*  os_swi3(8, 18, (int)draw_name, 0xAFF); */
  504.  
  505.   return(0);
  506. }
  507.  
  508.