6

I'm trying to make sure I understand my (digital) signal processing knowledge, by realizing a time-discrete version of a 1st order RC filter. (The background is that I'm trying to implement a PLL in software for SDR purposes, but this is a different story...)

My problem is that I thought I understood how to create the difference equation for such a filter, and therefore derive its coefficients. However when I plot the response in MATLAB using the freqz function - with the calculated a and b coefficients - I don't get what looks like an RC filter response.

I referenced the Wikipedia page on this topic (at http://en.wikipedia.org/wiki/Low-pass_filter#Discrete-time_realization), just to make sure I wasn't totally off in the weeds, but it still doesn't help. This details the difference equation as:

yi = alpha * xi + ( 1 - alpha ) * yi-1
where: alpha = sample period / ( RC + sample period )

An example:

fs = 96000.0;                         % Sample rate.
delta_t = 1.0 / fs;                   % Sample period.
fc = 5000.0;                          % Filter cut off frequency.
tau = 1 / ( 2 * pi * fc );            % Time constant of filter.
alpha = delta_t / ( tau + delta_t );  % Smoothing factor per Wikipedia page.
b = [ alpha ];                        % 'b' coefficients
a = [ 1.0, ( 1 - alpha ) ];           % 'a' coefficents
freqz( b, a, 1024, fs );              % 1024 point FFT used.

The result: enter image description here

Any thoughts as to where I'm going wrong? Have I totally misunderstood something?

Thanks in advance.

bitcyber
  • 131
  • 3
  • 11
  • In general, a first-order RC filter is a poor choice for DSP. The advantage when using discrete components is obvious -- fewer components. With DSP, that doesn't matter. – Ben Voigt Aug 04 '12 at 22:48
  • Yes - a good comment, and I couldn't agree with you more. It's just that much of the published PLL material assumes an analog implementation. There are of course many time discrete / software PLL articles published, but I want to make sure I understand the topic of PLLs totally. Thus I'm starting at the basics and came unstuck a little earlier than expected. – bitcyber Aug 04 '12 at 22:57

2 Answers2

3

You want your a(2) coefficient to be negative, since a represents the coefficient that appears on the left-hand side of the equation.

a(1)*y(n) + a(2)*y(n-1) - ... + a(na+1)*y(n-na) = b(1)*x(n) + b(2)*x(n-1) + ... + b(nb+1)*x(n-nb)

or equivalently,

a = a ./ a(1)
y(n) = b(1)*x(n) + b(2)*x(n-1) + ... + b(nb+1)*x(n-nb)
                 - a(2)*y(n-1) - ... - a(na+1)*y(n-na)

See the documentation for filter

With this correction, the response becomes

enter image description here

Ben Voigt
  • 277,958
  • 43
  • 419
  • 720
  • Thanks for this and an obvious oversight on my behalf. Much appreciated! I did try this, as I was playing around with the a coefficient values, but what led me to post my question was the -3 dB point of the resultant filter. It is quite a bit below 5000 Hz, around 4406 Hz for the -3.0103 dB point (i.e., 10 * log10[0.5]). Anyway, the main problem is solved. Thanks again. – bitcyber Aug 04 '12 at 22:57
2

Ben Voigt gave the correct answer. He just didn't show how to change bitcyber's code to make the code produce his graphs. I needed this StackExchange answer to make sense of what Ben said, and then I changed the third last line to "a = [ 1.0, ( alpha -1 ) ];". 

fs = 96000.0;                         % Sample rate.
delta_t = 1.0 / fs;                   % Sample period.
fc = 5000.0;                          % Filter cut off frequency.
tau = 1 / ( 2 * pi * fc );            % Time constant of filter.
alpha = delta_t / ( tau + delta_t );  % Smoothing factor per Wikipedia page.
b = [ alpha ];                        % 'b' coefficients
a = [ 1.0, ( alpha -1 ) ];           % 'a' coefficents
a = a ./ a(1);
freqz( b, a, 1024, fs );              % 1024 point FFT used.

Thanks for your answer Ben. It's just that I am a bit slower to catch on than most.

Community
  • 1
  • 1
MechtEngineer
  • 547
  • 6
  • 12