Temporary discrete implementation of the 1st order RC filter

Question

I try to make sure that I understand my (digital) knowledge of signal processing by implementing a discrete version of the 1st order RC filter. (It is assumed that I am trying to implement PLL in software for SDR purposes, but this is a completely different story ...)

My problem is that I thought I understood how to create a difference equation for such a filter, and therefore get its coefficients. However, when I draw the answer in MATLAB using the freqz function - with the calculated coefficients a and b , I do not get what looks like the answer of the RC filter.

I referred to the Wikipedia page on this topic ( http://en.wikipedia.org/wiki/Low-pass_filter#Discrete-time_realization ), just to make sure that I’m not completely at the weed, but that still doesn’t help. This describes the difference equation as:

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

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. 

Result:

Any thoughts on where I'm wrong? I didn’t understand something at all?

Thanks in advance.

Answer 1

You want your coefficient a(2) be negative, since a represents the coefficient that appears on the left 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 documentation for filter

---

With this correction, the answer becomes

Answer 2

Ben Voigt gave the correct answer. He simply did not show how to change the bitcyber code so that the code generated its graphics. I need this StackExchange answer to understand 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 reply Ben. It's just that I'm a little slower than most.