#if !JAVA /* * math.cs * Copyright (c) 2008-2009 kbinani * * This file is part of bocoree. * * bocoree is free software; you can redistribute it and/or * modify it under the terms of the BSD License. * * bocoree is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. */ using System; using System.Collections.Generic; namespace bocoree { public partial class math { private const double _PI2 = 2.0 * Math.PI; private const double _PI4 = 4.0 * Math.PI; private const double _PI6 = 6.0 * Math.PI; private const double _PI8 = 8.0 * Math.PI; public enum WindowFunctionType { Hamming, rectangular, Gauss, Hann, Blackman, Bartlett, Nuttall, Blackman_Harris, Blackman_Nattall, flap_top, Parzen, Akaike, Welch, Kaiser, } public static double window_func( WindowFunctionType type, double x ) { switch ( type ) { case WindowFunctionType.Akaike: return wnd_akaike( x ); case WindowFunctionType.Bartlett: return wnd_bartlett( x ); case WindowFunctionType.Blackman: return wnd_blackman( x ); case WindowFunctionType.Blackman_Harris: return wnd_blackman_harris( x ); case WindowFunctionType.Blackman_Nattall: return wnd_blackman_nattall( x ); case WindowFunctionType.flap_top: return wnd_flap_top( x ); case WindowFunctionType.Gauss: throw new ApplicationException( "too few argument for Gauss window function" ); case WindowFunctionType.Hamming: return wnd_hamming( x ); case WindowFunctionType.Hann: return wnd_hann( x ); case WindowFunctionType.Kaiser: throw new ApplicationException( "too few argument for Kaiser window function" ); case WindowFunctionType.Nuttall: return wnd_nuttall( x ); case WindowFunctionType.Parzen: return wnd_parzen( x ); case WindowFunctionType.rectangular: return wnd_rectangular( x ); case WindowFunctionType.Welch: return wnd_welch( x ); } return 0.0; } public static double window_func( WindowFunctionType type, double x, params double[] param ) { switch ( type ) { case WindowFunctionType.Akaike: return wnd_akaike( x ); case WindowFunctionType.Bartlett: return wnd_bartlett( x ); case WindowFunctionType.Blackman: return wnd_blackman( x ); case WindowFunctionType.Blackman_Harris: return wnd_blackman_harris( x ); case WindowFunctionType.Blackman_Nattall: return wnd_blackman_nattall( x ); case WindowFunctionType.flap_top: return wnd_flap_top( x ); case WindowFunctionType.Gauss: return wnd_gauss( x, param[0] ); case WindowFunctionType.Hamming: return wnd_hamming( x ); case WindowFunctionType.Hann: return wnd_hann( x ); case WindowFunctionType.Kaiser: return wnd_kaiser( x, param[0] ); case WindowFunctionType.Nuttall: return wnd_nuttall( x ); case WindowFunctionType.Parzen: return wnd_parzen( x ); case WindowFunctionType.rectangular: return wnd_rectangular( x ); case WindowFunctionType.Welch: return wnd_welch( x ); } return 0.0; } public static double wnd_kaiser( double x, double alpha ) { if ( 0.0 <= x && x <= 1.0 ) { double t = 2.0 * x - 1.0; return besi0( Math.PI * alpha * Math.Sqrt( 1.0 - t * t ) ) / besi0( Math.PI * alpha ); } else { return 0.0; } } public static double wnd_welch( double x ) { if ( 0.0 <= x && x <= 1.0 ) { return 4.0 * x * (1.0 - x); } else { return 0.0; } } public static double wnd_akaike( double x ) { if ( 0.0 <= x && x <= 1.0 ) { return 0.625 - 0.5 * Math.Cos( _PI2 * x ) - 0.125 * Math.Cos( _PI4 * x ); } else { return 0.0; } } public static double wnd_parzen( double x ) { double x0 = Math.Abs( x ); if ( x0 <= 1.0 ) { return (0.75 * x0 - 1.5) * x0 * x0 + 1.0; } else { x0 = 2.0 - x0; return 0.25 * x0 * x0 * x0; } } public static double wnd_flap_top( double x ) { if ( 0.0 <= x && x <= 1.0 ) { return 1.0 - 1.93 * Math.Cos( _PI2 * x ) + 1.29 * Math.Cos( _PI4 * x ) - 0.388 * Math.Cos( _PI6 * x ) + 0.032 * Math.Cos( _PI8 * x ); } else { return 0.0; } } public static double wnd_blackman_nattall( double x ) { if ( 0.0 <= x && x <= 1.0 ) { return 0.3635819 - 0.4891775 * Math.Cos( _PI2 * x ) + 0.1365995 * Math.Cos( _PI4 * x ) - 0.0106411 * Math.Cos( _PI6 * x ); } else { return 0.0; } } public static double wnd_blackman_harris( double x ) { if ( 0.0 <= x && x <= 1.0 ) { return 0.35875 - 0.48829 * Math.Cos( _PI2 * x ) + 0.14128 * Math.Cos( _PI4 * x ) - 0.01168 * Math.Cos( _PI6 * x ); } else { return 0.0; } } public static double wnd_nuttall( double x ) { if ( 0.0 <= x && x <= 1.0 ) { return 0.355768 - 0.487396 * Math.Cos( _PI2 * x ) + 0.144232 * Math.Cos( _PI4 * x ) - 0.012604 * Math.Cos( _PI6 * x ); } else { return 0.0; } } public static double wnd_bartlett( double x ) { if ( 0.0 <= x && x <= 1.0 ) { return 1.0 - 2.0 * Math.Abs( x - 0.5 ); } else { return 0.0; } } public static double wnd_blackman( double x ) { if ( 0.0 <= x && x <= 1.0 ) { return 0.42 - 0.5 * Math.Cos( _PI2 * x ) + 0.08 * Math.Cos( _PI4 * x ); } else { return 0.0; } } public static double wnd_hann( double x ) { if ( 0.0 <= x && x <= 1.0 ) { return 0.5 - 0.5 * Math.Cos( _PI2 * x ); } else { return 0.0; } } public static double wnd_gauss( double x, double sigma ) { return Math.Exp( -x * x / (sigma * sigma) ); } public static double wnd_rectangular( double x ) { if ( 0.0 <= x && x <= 1.0 ) { return 1.0; } else { return 0.0; } } public static double wnd_hamming( double x ) { if ( 0.0 <= x && x <= 1.0 ) { return 0.54 - 0.46 * Math.Cos( _PI2 * x ); } else { return 0.0; } } public static double besi0( double x ) { int i; double w, wx375; double[] a = { 1.0, 3.5156229, 3.0899424, 1.2067492, 0.2659732, 0.0360768}; double[] b = { 0.39894228, 0.013285917, 0.002253187, -0.001575649, 0.009162808, -0.020577063, 0.026355372, -0.016476329}; if ( x < 0.0 ) { return 0.0; } if ( x <= 3.75 ) { wx375 = x * x / 14.0625; w = 0.0045813; for ( i = 5; i >= 0; i-- ) { w = w * wx375 + a[i]; } return w; } wx375 = 3.75 / x; w = 0.003923767; for ( i = 7; i >= 0; i-- ) { w = w * wx375 + b[i]; } return w / Math.Sqrt( x ) * Math.Exp( x ); } public static double erfcc( double x ) { double t, z, res; z = Math.Abs( x ); t = 1.0 / (1.0 + 0.5 * z); res = t * Math.Exp( -z * z - 1.26551223 + t * (1.00002368 + t * (0.37409196 + t * (0.09678418 + t * (-0.18628806 + t * (0.27886807 + t * (-1.13520398 + t * (1.48851587 + t * (-0.82215223 + t * 0.17087277)))))))) ); if ( x < 0.0 ) { res = 2.0 - res; } return res; } } } #endif