I've two signals X and Y. These are proximity sensors reading measuring vibrations of a rotating shaft. Now I want to plot orbit plots for these signals. I can plot unfiltered orbits with just X and Y. But since there is lot of noise in the signal I filtered it. Now there are dominant 1X frequency component and may be other sub frequency components. I want to plot the filtered orbit of 1X say 656 Hz. The sampling frequency I used is 1000000 Hz. I know band pass filter could give me only the 1X component. But I don't have DSP toolbox. Could someone help me how to design and implement narrow band pass filter without toolbox. Any suggestion or hints or idea would be helpful. Thanks in advance.

  • $\begingroup$ Welcome to DSP.SE! I'd suggest getting yourself a toolbox: Scilab, Octave, or Python are popular (and free) tools. $\endgroup$ – Peter K. Mar 16 '15 at 12:17
  • $\begingroup$ I have all of them in my Personal PC. But this is part of my Project at Univ. So I need to implement in MATLAB only which unfortunately doesn't have DSP toolbox. $\endgroup$ – Agni Mar 16 '15 at 12:34
  • $\begingroup$ OK! Thanks for the clarification. Let's see if we can come up with a good answer. Are you sure the sampling frequency is 1MHz? that seems like overkill if your signal of interest is only 656Hz. $\endgroup$ – Peter K. Mar 16 '15 at 12:36
  • $\begingroup$ I've data set coming from another client. They mentioned in their data sheets the sampling rate as 2MHz and since the shaft is rotating at 38 K rpm( 633.33 Hz). 656 Hz is the fundamental frequency component (1X). I tried decreasing it earlier but I#m loosing resolution and the plots are going weird. The total data set has 40040 data points $\endgroup$ – Agni Mar 16 '15 at 12:47
  • $\begingroup$ Ah. Just because your shaft speed is 38,000 RPM, does not mean the frequency of interest is necessarily 633.33Hz. It will depend on the shaft geometry and bearing types. Do you know what the frequency range of the noise you want to get rid of is? I'd suggest designing the filter in python (or scilab) and just transferring the coefficients to matlab and using filter to apply them. $\endgroup$ – Peter K. Mar 16 '15 at 20:15

The steps that I would follow to handcraft a bandpass filter are these:

  1. Figure out how wide (in Hz) you want the filter to be.
  2. Create a low-pass filter with half the one-sided bandwidth as your desired filter bandwidth. You can calculate the sinc function with basic functions, since $sinc(x) = \frac{sin(\pi x)}{\pi x}$. You can adjust the bandwidth by adjusting the time scale in your sinc function. If n is the set of all integers between -N and N, inclusive, and you make your filter taps equal to $sinc(bn)$, then the filter's normalized bandwidth will be equal to $b$.
  3. Modulate the lowpass filter to the center frequency you want for your bandpass filter. You do that by multiplying it, element by element, with a sine wave of that frequency.
  • $\begingroup$ Thanks for your answer. If I had a DSP toolbox how would I implement it? I tried it with basic filters like band pass ( butterworth) but it seem to fail in that context too. Could you please look into this. stackoverflow.com/questions/29111768/… $\endgroup$ – Agni Mar 19 '15 at 9:09

Assuming 2 MHz sampling rate, by just simple addition of every 20 samples you get a low pass filter that will keep the noise down and will still allow you to recover every problem related to the 38K RPM source signal which you could cope with successfully since 100K is maintaining the Nyquist requirements.


If you'd like to design IIR, use pole-zero placement method. That's what people used before dsp toolboxes emerged.


A crude approach do design an IIR filter for a single frequency $\omega, 0 \le \omega \le 1$ would be as follows:

  • Placing poles at $ \exp(\pm j \pi \omega)$ would give infinite gain for that frequency, so move them "a bit" inside to a radius $< 1$, say, $r$.

  • Place zeros at $\pm 1$ to "guarantee" a bandpass filter

  • Place zeros at other positions on the unit circle to increase the attenuation

  • Figure out the gain at $\omega$, where it would/should be the highest, and compensate for that

For your example:

omega = 656/(1000000)/2;
r = 0.99;
b = conv([1 -1], [1 1]);
a = conv([1 -r*exp(j*pi*omega)],[1 -r*exp(-j*pi*omega)]);

will form a start. Considering that the 1X frequency is rather low, you could add a number of zeros at some higher frequencies:

zerof = 0.1; % Corresponding to 50000 Hz
b = conv(b, conv([1 -exp(j*pi*zerof)],[1 -exp(-j*pi*zerof)]));
zerof = 0.5; % Corresponding to 250000 Hz

and multiple poles at the wanted frequency;

a = conv(a,a);

Once you are happy, you just need to figure out the gain and compensate for that.

This filter has about 80 dB attenuation for frequencies of 50000 Hz and higher (as someone pointed out, the sampling frequency is really high if correct).

  • $\begingroup$ This would be a quite rudimentary and manual pole-zero placement approach (there are automatic methods which are able to find better positions than those I more or less just made up). $\endgroup$ – Oscar Mar 19 '15 at 19:32

If you know how to design and apply Kaiser window, i would aim for a brickwall bandpass that might have (symmetrical and non-causal) impulse response

$$ h[n] = \cos(\omega_0 n) \cdot \frac{\sin(\omega_1 n/2)}{\omega_1 n/2} \quad \forall n \in \mathbb{Z} $$

assuming $0 < \omega_1 < \omega_0$, the bandedges should be at $\omega_0 \pm \omega_1$. the number of samples of $h[n]$ to keep and the "$\beta$ factor" for the Kaiser window would be things you have to fiddle with.


Your Answer

