Time discrete implementation of 1st order RC filter

Question

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:

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

Thanks in advance.

Answer 1

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

Answer 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.