timidity/filter.c
author Sam Lantinga <slouken@libsdl.org>
Sat, 01 Jun 2013 22:30:04 -0700
changeset 627 438fb12df34f
parent 518 8bc9b5fd2aae
permissions -rw-r--r--
Updated iOS project to include OGG support directly
     1 /*
     2     TiMidity -- Experimental MIDI to WAVE converter
     3     Copyright (C) 1995 Tuukka Toivonen <toivonen@clinet.fi>
     4 
     5     This program is free software; you can redistribute it and/or modify
     6     it under the terms of the Perl Artistic License, available in COPYING.
     7 
     8    filter.c: written by Vincent Pagel ( pagel@loria.fr ) 
     9  
    10    implements fir antialiasing filter : should help when setting sample
    11    rates as low as 8Khz.
    12 
    13    April 95
    14       - first draft
    15 
    16    22/5/95
    17       - modify "filter" so that it simulate leading and trailing 0 in the buffer
    18    */
    19 
    20 #include <stdio.h>
    21 #include <string.h>
    22 #include <math.h>
    23 #include <stdlib.h>
    24 #include "config.h"
    25 #include "common.h"
    26 #include "ctrlmode.h"
    27 #include "instrum.h"
    28 #include "filter.h"
    29 
    30 /*  bessel  function   */
    31 static float ino(float x)
    32 {
    33     float y, de, e, sde;
    34     int i;
    35     
    36     y = x / 2;
    37     e = 1.0;
    38     de = 1.0;
    39     i = 1;
    40     do {
    41 	de = de * y / (float) i;
    42 	sde = de * de;
    43 	e += sde;
    44     } while (!( (e * 1.0e-08 - sde > 0) || (i++ > 25) ));
    45     return(e);
    46 }	
    47 
    48 /* Kaiser Window (symetric) */
    49 static void kaiser(float *w,int n,float beta)
    50 {
    51     float xind, xi;
    52     int i;
    53     
    54     xind = (float)((2*n - 1) * (2*n - 1));
    55     for (i =0; i<n ; i++) 
    56 	{
    57 	    xi = (float)(i + 0.5);
    58 	    w[i] = ino((float)(beta * sqrt((double)(1. - 4 * xi * xi / xind))))
    59 		/ ino((float)beta);
    60 	}
    61 }
    62 
    63 /*
    64  * fir coef in g, cuttoff frequency in fc
    65  */
    66 static void designfir(float *g , float fc)
    67 {
    68     int i;
    69     float xi, omega, att, beta ;
    70     float w[ORDER2];
    71     
    72     for (i =0; i < ORDER2 ;i++) 
    73 	{
    74 	    xi = (float) (i + 0.5);
    75 	    omega = (float)(PI * xi);
    76 	    g[i] = (float)(sin( (double) omega * fc) / omega);
    77 	}
    78     
    79     att = 40.; /* attenuation  in  db */
    80     beta = (float) (exp(log((double)0.58417 * (att - 20.96)) * 0.4) + 0.07886 
    81 	* (att - 20.96));
    82     kaiser( w, ORDER2, beta);
    83     
    84     /* Matrix product */
    85     for (i =0; i < ORDER2 ; i++)
    86 	g[i] = g[i] * w[i];
    87 }
    88 
    89 /*
    90  * FIR filtering -> apply the filter given by coef[] to the data buffer
    91  * Note that we simulate leading and trailing 0 at the border of the 
    92  * data buffer
    93  */
    94 static void filter(sample_t *result,sample_t *data, int32 length,float coef[])
    95 {
    96     int32 sample,i,sample_window;
    97     int16 peak = 0;
    98     float sum;
    99 
   100     /* Simulate leading 0 at the begining of the buffer */
   101      for (sample = 0; sample < ORDER2 ; sample++ )
   102 	{
   103 	    sum = 0.0;
   104 	    sample_window= sample - ORDER2;
   105 	   
   106 	    for (i = 0; i < ORDER ;i++) 
   107 		sum += (float)(coef[i] *
   108 		    ((sample_window<0)? 0.0 : data[sample_window++])) ;
   109 	    
   110 	    /* Saturation ??? */
   111 	    if (sum> 32767.) { sum=32767.; peak++; }
   112 	    if (sum< -32768.) { sum=-32768; peak++; }
   113 	    result[sample] = (sample_t) sum;
   114 	}
   115 
   116     /* The core of the buffer  */
   117     for (sample = ORDER2; sample < length - ORDER + ORDER2 ; sample++ )
   118 	{
   119 	    sum = 0.0;
   120 	    sample_window= sample - ORDER2;
   121 
   122 	    for (i = 0; i < ORDER ;i++) 
   123 		sum += data[sample_window++] * coef[i];
   124 	    
   125 	    /* Saturation ??? */
   126 	    if (sum> 32767.) { sum=32767.; peak++; }
   127 	    if (sum< -32768.) { sum=-32768; peak++; }
   128 	    result[sample] = (sample_t) sum;
   129 	}
   130     
   131     /* Simulate 0 at the end of the buffer */
   132     for (sample = length - ORDER + ORDER2; sample < length ; sample++ )
   133 	{
   134 	    sum = 0.0;
   135 	    sample_window= sample - ORDER2;
   136 	    
   137 	    for (i = 0; i < ORDER ;i++) 
   138 		sum += (float)(coef[i] *
   139 		    ((sample_window>=length)? 0.0 : data[sample_window++])) ;
   140 	    
   141 	    /* Saturation ??? */
   142 	    if (sum> 32767.) { sum=32767.; peak++; }
   143 	    if (sum< -32768.) { sum=-32768; peak++; }
   144 	    result[sample] = (sample_t) sum;
   145 	}
   146 
   147     if (peak)
   148 	ctl->cmsg(CMSG_ERROR, VERB_NORMAL, 
   149 		  "Saturation %2.3f %%.", 100.0*peak/ (float) length);
   150 }
   151 
   152 /***********************************************************************/
   153 /* Prevent aliasing by filtering any freq above the output_rate        */
   154 /*                                                                     */
   155 /* I don't worry about looping point -> they will remain soft if they  */
   156 /* were already                                                        */
   157 /***********************************************************************/
   158 void antialiasing(Sample *sp, int32 output_rate )
   159 {
   160     sample_t *temp;
   161     int i;
   162     float fir_symetric[ORDER];
   163     float fir_coef[ORDER2];
   164     float freq_cut;  /* cutoff frequency [0..1.0] FREQ_CUT/SAMP_FREQ*/
   165  
   166 
   167     ctl->cmsg(CMSG_INFO, VERB_NOISY, "Antialiasing: Fsample=%iKHz",
   168 	      sp->sample_rate);
   169  
   170     /* No oversampling  */
   171     if (output_rate>=sp->sample_rate)
   172 	return;
   173     
   174     freq_cut= (float) output_rate / (float) sp->sample_rate;
   175     ctl->cmsg(CMSG_INFO, VERB_NOISY, "Antialiasing: cutoff=%f%%",
   176 	      freq_cut*100.);
   177 
   178     designfir(fir_coef,freq_cut);
   179     
   180     /* Make the filter symetric */
   181     for (i = 0 ; i<ORDER2 ;i++) 
   182 	fir_symetric[ORDER-1 - i] = fir_symetric[i] = fir_coef[ORDER2-1 - i];
   183     
   184     /* We apply the filter we have designed on a copy of the patch */
   185     temp = safe_malloc(sp->data_length);
   186     memcpy(temp,sp->data,sp->data_length);
   187     
   188     filter(sp->data,temp,sp->data_length/sizeof(sample_t),fir_symetric);
   189     
   190     free(temp);
   191 }