1
$\begingroup$

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

I followed these two links http://doctord.webhop.net/Courses/textbooks/Smith_DSP/dsp_book_Ch16.pdf and https://tomroelandts.com/articles/how-to-create-a-simple-low-pass-filter

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): enter image description here Filtering to 4th harmonic: enter image description here Filtering to 10kh: enter image description here

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 license header, choose License Headers in Project Properties.
 * 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 license header, choose License Headers in Project Properties.
 * 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);
            AdsrEnvelope ae = new AdsrEnvelope(1, 0.05, 0.2, true, 0.8, 0.25, sampleRate);
            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 = ae.getADSRvolume(exportIndex);
                    //pointVolume *= t.drawTremulantPoint(pointVolume, ae.getSustainStartPoint());
                    sg.setAmplitude(pointVolume);
                    points[exportIndex] = point;
                    points[exportIndex] = al.calculateAmpLimit(points[exportIndex]);
                } else {
                    exportIndex--;
                }
                /*pointVolume = ae.getADSRvolume(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);
        }

    }

}

$\endgroup$
6
  • 1
    $\begingroup$ 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. $\endgroup$ – Bob Mar 5 at 8:39
  • $\begingroup$ 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. $\endgroup$ – Edward Eddy67716 Mar 5 at 21:27
  • $\begingroup$ 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. $\endgroup$ – a concerned citizen Mar 6 at 8:08
  • $\begingroup$ Could you please explain using the variables I have in my code? $\endgroup$ – Edward Eddy67716 Mar 8 at 6:41
  • $\begingroup$ @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? $\endgroup$ – Edward Eddy67716 Mar 9 at 19:53
0
$\begingroup$

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)));

test

$\endgroup$
9
  • $\begingroup$ 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. $\endgroup$ – Edward Eddy67716 Mar 10 at 21:15
  • $\begingroup$ @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. $\endgroup$ – a concerned citizen Mar 10 at 21:31
  • $\begingroup$ Actually it is not that noticeable. $\endgroup$ – Edward Eddy67716 Mar 10 at 21:36
  • $\begingroup$ Is there a way I can apply resonance to this filter like one can to a Two pole filter? $\endgroup$ – Edward Eddy67716 Mar 11 at 6:32
  • $\begingroup$ @EdwardEddy67716 Sounds like you need an IIR. $\endgroup$ – a concerned citizen Mar 11 at 8:14
0
$\begingroup$

I have managed to get the whole thing to work.

Code to do so:

/**
 *
 * @author Edward Jenkins
 * @version 1.2
 */
public abstract class SincSpecs implements ISimpleFilter{
    
    // instance variables
    private double cutoffFrequency;
    private double sampleRate;
    private double resonance;
    private boolean highPass;
    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;
    //private int index;
    
    // constructor
    public SincSpecs(double cutoffFrequency, double sampleRate, byte resonance,
            boolean highPass, int band, ISample sample) {
        this.cutoffFrequency = cutoffFrequency;
        this.resonance = (double)resonance / 127;
        this.highPass = highPass;
        this.sampleRate = sampleRate;
        cutoffAmount = cutoffFrequency  / this.sampleRate;
        transitionBand = band / sampleRate;
        //index = 0;
        n = (int)(Math.ceil(4 / transitionBand));
        
        // make sure length is odd
        if (n % 2 == 0) {
            n += 1;
        }
        order = n;
        midPoint = order / 2 + 1;
        impulseResponce = new double[n];
        sumnatedValue = 0;

        updateImpulseResponce();

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

        double inputPoint = sample.generateWaveformPoint();

        while (!canOutput()) {
            inputFilterPoint(inputPoint);
        }
    }
    
    public double getTransitionBand() {
        return transitionBand;
    }
    
    @Override
    public void setCutoff(double cutoffFrequency) {
        
        this.cutoffFrequency = cutoffFrequency;
        cutoffAmount = this.cutoffFrequency / sampleRate;
        updateImpulseResponce();
    }
    
    public void updateImpulseResponce() {
        
        // 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;
        }
        
        // invert the impulse responce if high pass is true
        if (this.highPass) {
            invertImpulseForHighPass();
        }
    }
    
    // input points
    protected void inputFilterPoint(double point) {

        sampleCache[sampleCacheLength] = point;
        sampleCacheLength++;

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

    protected void incrementCache(double value) {

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

    // can output
    protected boolean canOutput() {
        boolean result = false;

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

        return result;
    }

    // sinc
    protected double sinc(double value) {
        if (value == 0) {
            value = 1;
        } else {
            value = Math.sin(value) / (value) + resonance *  Math.sin(value);
        }
        return value;
    }

    protected void invertImpulseForHighPass() {

        for (int i = 0; i < impulseResponce.length; i++) {
            impulseResponce[i] *= -1;
        }
        impulseResponce[midPoint] += 1;
    }

    @Override
    public double filter(double point) {

        //index++;
        
        value = 0;

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

        return value;
    }
}
```
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.