Creating a 50hz comb filter in MATLAB

Question

I have a signal in a 1d vector X, and I need to remove 50Hz and its resonant frequencies. I've looked into using fdesign.comb but I'm not sure how to apply that filter on my data (after building it using design).

Alternatively, I think I can also do this using iircomb but again, I'm not sure how to exactly. If someone could step me through this (and how to then apply it to a signal X) it would be appreciated.

Answer 1

In order to use fdesign.comb you have to do something like the following to specify the parameters of the filter

D = fdesign.comb('notch','N,BW',8,.01); 

but then you must generate the filter as you mentioned in your question.

H = design(D); % H is a DFILT

you can then view the filter like this

fvtool(H)

One way to apply the filter to your data is to use the filter function

y=filter(H.Numerator,H.Denominator,inArray)

My favorite way to design a filter I've never used is with fdatool because you can play around with all of your parameters. Here is how I would do it:

• Type fdatool at the command prompt and a window pops up
• Specify Notching in the dropdown that initially says Differentiator
• Inder IIR choose comb
• Input your parameters
• click Design Filter

Once you play around and it looks as you wish, you can either export the coefficients or generate the code that would create this filter.

To export coefficients

• go to file and click export
• make note of the variable names and click export
• now you can use the filter command with y = filter(num,den,data);

To generate code

• go to file and click generate MATLAB code
• save the file

You see in this file something like the following

Fs = 48000;  % Sampling Frequency

N     = 10;    % Order
BW    = 1200;  % Bandwidth
Apass = 1;     % Bandwidth Attenuation

[b, a] = iircomb(N, BW/(Fs/2), Apass);
Hd     = dfilt.df2(b, a);

the line with iircomb is what generates the filter coefficients, you can now use y=filter(b,a,data);

Hope this helps.

Answer 2

I needed to do something like this just recently!

Here is some example Matlab code that you might like to use as a starting point . . .

clear all; 
%% // Make a harmonic signal and add it to another sine
sr = 44100; 
dt = 1/sr;
tAxis = dt:dt:2;

freq = 50;
x = mean(sin (2*pi*freq*(1:10)'*tAxis))'; %Harmonic noise
x = x+ sin(2*pi*440*tAxis)';


%% // Fractional delay comb filter
rho = 0.95; %This sets how narrow the rejection is. Must be between 0 and 1

delay = (sr/2)/freq - 1;
N = 2*floor(delay);



%// LAGRANGE  h=lagrange(N,delay) returns order N FIR 
%//          filter h which implements given delay 
%//          (in samples).  For best results, 
%//          delay should be near N/2 +/- 1.
n = 0:N;
h = ones(1,N+1);
for k = 0:N
    index = find(n ~= k);
    h(index) = h(index) *  (delay-k)./ (n(index)-k);
end


b = [1 -h];
a = [1 -rho*h];

nfft = 2^14;
figure; freqz(b, a, nfft, 'half', sr)

xf = filter(b,a,x);

soundsc([x; xf],sr);

The script makes a 440 Hz sinusoid and adds 50 Hz harmonic noise. A 50 Hz comb filter is created using a fractional delay filter. THe signal is then filtered with the filter and the output is played as a sound before and after the filtering process.

For more detail try looking here

http://ieeexplore.ieee.org/xpl/freeabs_all.jsp?arnumber=951433

Hope that this helps.