# How do I sort out my window sinc low pass filter

I started writing a windowed sinc low pass filter, but I cannot seam to get the code figured out right.

Thanks to @aconcerned It now only has one problem, The waveforms are not normalised and I would like them to be. I want them to maintain a somewhat constant amplitude no matter how much they filter.

Filtering to fundamental (I expect more falloff here due to band limit): Filtering to 4th harmonic: Filtering to 10kh:

I tried this filter on square and sawtooth waves so far. I am building it for resampling on a new music tracker format that I am coding.

Filter

/*
* To change this template file, choose Tools | Templates
* and open the template in the editor.
*/
package effects;

/**
*
* @author Edward Jenkins
*/
public class SincLowPassFilter {

public static final double DEF_BAND = 0;

// instance variables
private double cutoffFrequency;
private double sampleRate;
private double value;
private double sumnatedValue;
private double cutoffAmount;
private double transitionBand;
private double window;
private double[] impulseResponce;
private double[] sampleCache;
private int sampleCacheLength;
private int n;
private int order;
private int midPoint;

// constructor
public SincLowPassFilter(double cutoffFrequency, double sampleRate,
int band) {
this.cutoffFrequency = cutoffFrequency;
cutoffAmount = cutoffFrequency  / sampleRate;
transitionBand = band / sampleRate;
n = (int)(Math.ceil(4 / transitionBand));

// make sure length is odd
if (n % 2 == 0) {
n += 1;
}
order = n - 1;
midPoint = order / 2;
impulseResponce = new double[n];
sumnatedValue = 0;

// get window of filtering
for (int i = 0; i < n; i++) {

impulseResponce[i] = cutoffAmount
* sinc(2 * cutoffAmount * Math.PI * (i - midPoint));

impulseResponce[midPoint] = cutoffAmount;

window = 0.54 - 0.46 * Math.cos(2 * Math.PI * i / order);

impulseResponce[i] *= window;
}

// sumnate all filter kernal values
double sum = 0;

for (int i = 0; i < n; i++) {
sum += impulseResponce[i];
}

for (int i = 0; i < n; i++) {
impulseResponce[i] /= sum;
impulseResponce[i] *= 1.2;
}

sampleCache = new double[n];
sampleCacheLength = 0;
}

// low pass filter
public void inputPoint(double point) {

value = 0;
/*value = impulseResponce[sampleCacheLength];

if (sampleCacheLength >= n-1) {
sampleCacheLength = 0;
} else {
sampleCacheLength++;
}*/

if (sampleCacheLength < n) {

sampleCache[sampleCacheLength] = point;
sampleCacheLength++;
} else {
for (int i = 0, j = n - 1; i < n; i++, j--) {
value += sampleCache[j] * impulseResponce[i];
}
incrementCache(point);
}
}

private void incrementCache(double value) {

for(int i = 0; i < sampleCacheLength - 1; i++) {
sampleCache[i] = sampleCache[i + 1];
}
sampleCache[sampleCacheLength - 1] = value;
}

// can low pass
public boolean canLowPass() {
boolean result = false;

if (sampleCacheLength >= n) {
result = true;
}

return result;
}

// sinc
public double sinc(double value) {
if (value == 0) {
value = 1;
} else {
value = Math.sin(value) / value;
}
return value;
}

public double lowPass() {
return value;
}
}



Execution

/*
* To change this template file, choose Tools | Templates
* and open the template in the editor.
*/
package sound.generator;

import effects.*;
import java.io.IOException;
import sound.generator.waveforms.*;
import sound.formatter.*;
import music.note.NoteRange;

/**
*
* @author Edward Jenkins
*/
public class SoundGenerator {

/**
* @param args the command line arguments
*/
public static void main(String[] args) {

NoteRange noteRange = new NoteRange(451, "C-1", "C9");
double frequency = noteRange.getFrequency("C4");
int sampleRate = 44100;
int bitRate = 16;
byte[] outputBytes;
double duration = 5;

// harmoncis
double[] harmonicVolumes = {1, 0, 1, 0, 0.2, 0, 0.2, 0, 0.09, 0, 0.08, 0, 0.07};

SawtoothGenerator sg = new SawtoothGenerator(frequency, 0);
try {
Celeste c = new Celeste(1, true, true, true, 0.75, sg);
Vibrato v = new Vibrato(0.66, 0.125, frequency, 0, sampleRate);
Tremulant t = new Tremulant(0.66, 0.75, 0.5, sampleRate, false);
SincLowPassFilter lpf = new SincLowPassFilter(261.63 * 6, sampleRate, 400);
int sampleLength = (int) (duration * sampleRate);
AmpLimiter al = new AmpLimiter(bitRate, false);
double point;
double[] points = new double[sampleLength];
double pointVolume;

for (int exportIndex = 0; exportIndex < sampleLength; exportIndex++) {
if (exportIndex == sampleRate * 4) {
ae.setSustain(false);
}
//sg.setFrequency(v.getVibratoFrequency(exportIndex));
point = c.detune();
lpf.inputPoint(point);
if (lpf.canLowPass()) {
point = lpf.lowPass();
//pointVolume *= t.drawTremulantPoint(pointVolume, ae.getSustainStartPoint());
sg.setAmplitude(pointVolume);
points[exportIndex] = point;
points[exportIndex] = al.calculateAmpLimit(points[exportIndex]);
} else {
exportIndex--;
}
//pointVolume *= t.drawTremulantPoint(pointVolume, ae.getSustainStartPoint());
c.setAmplitudes(pointVolume);
points[exportIndex] = point;
points[exportIndex] = al.calculateAmpLimit(points[exportIndex]);*/

}

IWriteWave ww = new WriteWave(1, 44100, bitRate, points, "LPF Sawtooth C4",
false);

try {
ww.writeWav();
} catch (IOException e) {
e.printStackTrace();
}
} catch (Exception e) {
System.out.println(e);
}

}

}


• Despite my efforts I could not understand your problem. But I see one potential error, passing zero to sinc as far as I know (n - 1) / 2 will compute an integer division, and sinc(0) = 1. – Bob Mar 5 at 8:39
• I have deduced that the either the sinc or the window function are wrong. also if I set N to 100. I can get better results. – Edward Eddy67716 Mar 5 at 21:27
• The function generating the impulse response is wrong. After calculating it, try setting h[floor(N/2)] *= 2. Either that or halven the others, I can't tell which one is wrong, but that's where the difference is. Also, not sure if it's halvened, but you can easily check that h[floor(N/2)] should be equal to wp. If you have acces to the raw sinc() function, use it as wp*sinc(wp*k). Window optional. – a concerned citizen Mar 6 at 8:08
• Could you please explain using the variables I have in my code? – Edward Eddy67716 Mar 8 at 6:41
• @a concerned citizen Do you mean by my code, filterKernal[(int)Math.floor(n / 2)] *= 2;? Also should it be done before or after the sumnation? – Edward Eddy67716 Mar 9 at 19:53

I'm not good in Java, but it looks like you defined your sinc() below with sin(x)/x, no checks, then you are using it in your calculations as:

h[n]=sinc(2*f*n)=sin(2*f*n)/(2*f*n)


Then you are using i-(n-1)/2 for your sinc but the check is for i-n/2. And your window is calculated as:

w[n]=0.54-0.46*cos(2*pi*n*(N-1))


I won't question your way of calculating the order, but if I got it right and this is how you're calculating them, you're not doing a good job. What you call the filter "kernal" (it's spelled kernel, BTW) is the impulse response. This, and the window (Hamming, by the looks of it) are calculated as:

\begin{align} h[n]&=\omega_p\mathrm{sinc}(\pi\omega_p(n-M))=\dfrac{\sin(\omega_p\pi(n-M))}{\pi(n-M)}\tag{1} \\ w[n]&=0.54-0.46\cos(2\pi\dfrac{n}{N})\tag{2} \end{align}

Where $$\omega_p$$ is the normalized corner (cut-off) frequency, $$N$$ is the order, $$M=N/2$$, and $$n=0,\,1,\,...\,,\,N$$. With these you can define sinc as:

if(x == 0)
return 1.0;
else
return sin(pi*x)/(pi*x);


and use it as wp*sinc(wp*n) to calculate the impulse response. Then, the window:

0.54-0.46*cos(2*pi*x/N)


where N is the order (how did you come up with that formula for calculating the order?). Note, not the length (L, number of taps), not the middle point (M, order divided by 2): N=L-1, M=N/2.

To test this, I used Octave (% denote a comment):

L=33;         % length of the filter
N=L-1;        % order
M=N/2;        % mid point
n=[0:N];
w0=2;         % sampling frequency
wp=200/512;   % corner (cut-off) frequency, chosen to match the grid, below
% impulse response, h[n]
% h=wp*sinc(wp*(n-M));   % this is the built-in function
h=sin(wp*pi*(n-M))./(pi*(n-M));
h(M+1)=wp;   % NOTE: Octave starts the index at 1, for you it will be h[M]
% Hamming window
\$ w=hamming(L);   % this is the buit-in function
w=0.54-0.46*cos(2*pi*n/N);
w(M+1)=1;   % for you it will be w[M], see the NOTE above
% test the frequency response
H=fft(h.*w,1024)(1:512);
subplot(1,2,1); plot(h,"",w); subplot(1,2,2); plot(20*log10(abs(H)));


• You code has solved most of the problems, but the waveform is not normalised. Before it maintained a constant amplitude, but not is gets softer the more you filter. – Edward Eddy67716 Mar 10 at 21:15
• @EdwardEddy67716 What do you mean it gets softer? Plot your impulse response in Audacious and perform an FFT on it. It should show a correct lowpass filter like in my image. If the results are the same, then the problem is elsewhere. Maybe in the code with the actual filtering. – a concerned citizen Mar 10 at 21:31
• Actually it is not that noticeable. – Edward Eddy67716 Mar 10 at 21:36
• Is there a way I can apply resonance to this filter like one can to a Two pole filter? – Edward Eddy67716 Mar 11 at 6:32
• @EdwardEddy67716 Sounds like you need an IIR. – a concerned citizen Mar 11 at 8:14