fir source code

1 module fir;
2 
3 import dplug.core.nogc;
4 
5 import std.math;
6 
7 import std.stdio;
8 // /**
9 //  * A fast method of convolution of a long signal with a FIR filter
10 //  */
11 // void overlapAdd(float[] impulseResponse, ref float[] inputSignal)
12 // {
13 //     int m = impulseResponse.length;
14 //     int nx = 
15 // }
16 
17 /**
18  * IDCT algorithm need to extract the impulse response from sampled frequency plot
19  */
20 float[] inverseDiscreteCosineTransformReal(float[] frequencyPlot)
21 {
22     int N = cast(int)frequencyPlot.length * 2;
23 
24     float[] a = mallocSlice!float(N);
25 
26     for(int n = 0; n < (N / 2); ++n)
27     {
28         float x = 0;
29         for(int i = 1; i < (N / 2); ++i)
30         {
31             x += abs(frequencyPlot[i] * abs(cos( abs(2 * PI * i *((n - (cast(float)N - 1) / 2)) / cast(float)N))));
32         }
33         a[n] = (1 / cast(float)N ) * (frequencyPlot[0] + 2 * x);
34     }
35 
36     // for(int n = ())
37 
38     // foreach(index, )
39 
40     return a;
41 }
42 
43 // /**
44 //  * IDCT algorithm need to extract the impulse response from sampled frequency plot
45 //  */
46 // float[] inverseDiscreteCosineTransformReal(float[] frequencyPlot)
47 // {
48 //     int N = cast(int)frequencyPlot.length * 2;
49 
50 //     float[] a = mallocSlice!float(N);
51 
52 //     for(int k = 0; k < (N / 2); ++k)
53 //     {
54 //         float an = 0;
55 //         for(int n = 0; n < (N / 2); ++n)
56 //         {
57 //             an += frequencyPlot[n] * cos((PI * k * (2 * n + 1)) / (2 * N));
58 //         }
59 //         a[k] = 2 * an;
60 //     }
61 
62 //     // for(int n = ())
63 
64 //     // foreach(index, )
65 
66 //     return a;
67 // }
68 
69 
70 unittest
71 {
72     import std.stdio;
73     import std.conv;
74 
75     float[] frequencyPlot = [1.0, 1.0, 1.0, 0.001, 0.001, 0.001, 0.001, 0.001];
76     float[] expected = [0.04858366, 0.00364087, -0.05199205, -0.07047625, -0.02194221, 0.08695625, 0.21101949, 0.29421023, 0.29421023, 0.21101949, 0.08695625, -0.02194221, -0.07047625, -0.05199205, 0.00364087, 0.04858366]; 
77     // float[] frequencyPlot = [0.1, 2.1, 0.3, 4.2];
78     // float[] expected = [ 0.95238737, -0.80969772,  0.7286317,  -0.82132135,  -0.82132135,  0.7286317, -0.80969772, 0.95238737];
79     float[] actual = inverseDiscreteCosineTransformReal(frequencyPlot);
80 
81     assert(actual.length == expected.length, "Expected actual.length to be " ~ to!string(expected.length) ~ `. Instead got ` ~ to!string(actual.length));
82     foreach(index, e; expected)
83     {
84         if(e != actual[index])
85         {
86             writeln("Index: " ~ to!string(index));
87             assert(e == actual[index], "Error: expected " ~ to!string(e) ~ " to equal " ~ to!string(actual[index]));
88         }
89     }
90 }