Using an algorithm that gives the fundamental frequency across time, example https://www.mathworks.com/help/audio/ref/pitch.html (or some other one?)

how would I then use that, to create a harmonic mask of nHarmonics based on that fundamental frequency? So the mask will contain the f0 + nHarmonics harmonics.

The point of this mask will be to then use it as a filter for the audio mixture.

I want to do this using MATLAB.

edit: I forgot to mention that the f0 is not a constant value, it is a ranged f0 track from fmin to fmax, for example dominant f0 from 100Hz to 800Hz. So the f0 track contains different values, not just one value. I can understand a comb filter working for one value, but this is a bit more complex, as the f0 track has multiple different values, so now each one of those would have to have nHarmonics multiples above the f0 for me to create this harmonic mask.

  • $\begingroup$ Hey Dan- interesting problem, hopefully you are not here hoping someone will code it for you. It would be better if you showed your attempts at how you would do this, and what your underlying (lower level) signal processing questions would be. $\endgroup$ – Dan Boschen Oct 20 '18 at 3:19
  • $\begingroup$ I can code this myself, I am just trying to structure this harmonic mask in a way that I can use it as a filter. So basically, given the output from the pitch function, to get the harmonics, I would probably have to do some multiplications of the fundamental and then sum it all up to form this mask, this is my current method, but I am just trying to see if someone has a better way of doing this. Or maybe even point me to a method somewhere that someone already did similar to my problem. $\endgroup$ – Dan Oct 20 '18 at 4:02
  • $\begingroup$ is it that you know the fundamental frequency of a waveform and you want a comb filter that will allow the harmonics to come through and kill everything else? or a comb filter that will null out the harmonics? is that what you want? $\endgroup$ – robert bristow-johnson Oct 20 '18 at 6:09
  • $\begingroup$ nulling out the harmonics OR passing them and kill everything else is the type of filter I am trying to create. You mention a comb filter, this is interesting, how would that work? Keep in mind, I don't necessarily want to null all harmonics, only the ones related to the most dominant source (f0) from the pitch tracking step. $\endgroup$ – Dan Oct 20 '18 at 14:35
  • $\begingroup$ Dan, you need to mention my name @robertbristow-johnson so that I am notified. I totally forgot I left this question and never seen that it was answered. Yes a comb filter can be tuned to exactly your fundamental $f_0$ and notch it out as well as ever integer harmonic. $\endgroup$ – robert bristow-johnson Dec 20 '18 at 0:11

See this post which describes creating a notch filter using a delay and add:

Moving Average for Notch filtering

The filter used is generically given as

$$H(z) = 1+z^{-N}$$

This results in passing DC and placing zeros at the normalized radian frequency locations given by:

$$\omega_n = \frac{2\pi n}{N} + \frac{\pi}{N}$$

for all integers n and $\omega_n$ in the range of 0 to $2\pi$ to correspond to the primary frequency range of DC to the sampling rate.

This comes directly out of factoring the polynomial given by H(z), resulting in zeros on the complex unit circle at the angular frequencies given above. Because of the extra $\pi/N$ rotation on the unit circle, this filter works well for filtering odd harmonics based on prudent choice of N and therefore the sampling rate as an integer multiple of the notch location.

If instead a delay and subtraction is used:

$$H(z) = 1-z^{-N}$$

The result would be a zero (notch) at DC and and placing zeros at the normalized radian frequency locations given by:

$$\omega_n = \frac{2\pi n}{N}$$

Which comes from factoring $H(z) = 1-z^{-N}$ which results in the Nth roots of unity (zeros evenly spaced around the unit circle). This approach conveniently allows for placing a notch at every harmonic of a selected frequency. To introduce the more advanced approaches given below, it is useful at this point to visually plot the zeros in the z domain for this case, using the example $H(z) = 1- z^{-8}$, along with the associated frequency response:

Delay and subtract comb filter

Both of these approaches are referred to as "Comb Filters" and notice in this case how the notch occurs at a fundamental frequency and every integer multiple of that frequency. The requirement here is the sampling rate be an integer multiple of the notch frequency. The filter has linear phase which is important if phase distortion is of concern, but a potential drawback is the gradual loss in frequency components at other than the distinct harmonics.

An improved approach for the amplitude response would be to extend the concept of the second order notch filter which I have derived at this link given below.

Transfer function of second order notch filter

This will result in much flatter amplitude response over the pass frequencies, but the tight notches come at the cost of increased phase distortion in proximity of the notches which may be an issue depending on the application.

To extend this to multiple notches at integer multiples note that the background of the 2nd order notch filter is the simpler DC notch filter given by

$$ H(z) = \frac{1+a}{2}\frac{(z-1)}{(z-a)} $$

Which has a zero at $z=1$ and a pole placed inside the unit circle at $z = a$ in close proximity to that zero; the closer the placement the tighter the bandwidth but also the more precision (bit size) is required to achieve the desired performance and maintain stability. Notice below the comparison to a single zero as done in the Comb Filter implementations (with $N=1$) to an example case of the DC notch filter ($a=0.95$ for the example shown):

DC nulling filters

At this point you may see where I am going with this. We can combine the concept of the DC Nulling filter with the comb filter to achieve tight notches at integer harmonic locations by placing poles in close proximity to the zeros given in the Comb Filter example as depicted in the graphic below.

Comb with tight notches

Following similar polynomial factoring relationships into the complex plane, this is realized with the following transfer function below:

$$H(z) = K\frac{1-z^{-N}}{1-a^Nz^{-N}}$$

The amplitude normalization factor K if desired is determined by solving for the value in the passpand $H(z=\omega_n)$ where $\omega_n = \frac{2\pi n}{N} + \frac{\pi}{N}$ and inverting this, which ends up being:

$$K = \frac{1+a^8}{2}$$

The result for N= 8 and a = 0.95 is shown in the plot below.

Comb with tight notches freq response

The possible challenges with this approach is the significant passband phase distortion as N increases, and the amount of precision required in implementation.

If the phase distortion was an issue and the notches need not be as tight as the previous IIR implementation can provide, then a third option is to combine the linear phase notch response given by $H(z) = 1-z^{-N}$ cascaded with a simple 3-tap linear phase compensator similar to what I detailed in this link (the compensator shown is actually a raised minus cosine response; a raised positive cosine response would be the approach for this case, following the same process outlined at that post.)

how to make CIC compensation filter

In particular see the plot in that post for the interpolated CIC compensator which is implmented by simple zero insert (meaning to interpolate by 8, simply insert 7 zeros between each filter coefficient in the 3 tap compensator). This won't achieve the flatness and sharp nulls as done with the IIR approach above but will have significantly improved flatness over a wider passband range than the simple notch filter approach AND with linear phase. The cascade of the simple delay and subtract comb filter with the interpolated 3 tap FIR compensator (given by the two responses shown below) results in a significantly flatter passband response with minimal added complexity.

DC nulling with raised cosine compensator

  • $\begingroup$ DanB, perhaps we need to show the other Dan how to make a comb filter (that's what i would call it, a "notch filter" is usually one with a notch at a single frequency with two zeros on the unit circle) tuned to any fundamental frequency. Even such that the period is not an integer number of samples. $\endgroup$ – robert bristow-johnson Dec 20 '18 at 0:16
  • $\begingroup$ I was thinking the same thing! Yes you should link him to this, I think I know the post you are referring to - just got too distracted. (Would really like to go to that jazz club in Waltham sometime however... life goes by...) $\endgroup$ – Dan Boschen Dec 20 '18 at 0:18
  • $\begingroup$ (Ah you are referring to this same post....thought you meant another recent one that wanted notches at harmonic frequency locations). But yes, get your point an optimized generic comb would be interesting. My first thought is arbitrary resampling of a comb structure such as this but I assume some interesting reductions could take place in an optimized approach. $\endgroup$ – Dan Boschen Dec 20 '18 at 0:22
  • $\begingroup$ @robertbristow-johnson I guess this post almost covers the arbitrary frequency locations: dsp.stackexchange.com/questions/52357/… So we just implement the time delay using a fractional delay filter; where any arbitrary delay can be implemented. $\endgroup$ – Dan Boschen Dec 20 '18 at 0:41
  • 1
    $\begingroup$ to tune a comb to any arbitrary pitch, you need to make a good fractional delay element. to make an "optimized" generic comb, we need to consider "simply" designing a digital filter (IIR or FIR) and replacing the $z^{-1}$ elements with $e^{-s\tau}$ elements where $\tau$ is a precision delay (not limited to an integer multiple of the sampling period). you can design the teeth of your comb filter to have a Butterworth response (very flat teeth). $\endgroup$ – robert bristow-johnson Dec 20 '18 at 0:47

Your Answer

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

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