matlab filter() function artifacts

Hi all,

I am facing the following problem:

a 100 Hz sine tone is fed from the input to the output of a system. When I am directly passing the input to the output and check the spectrum then the input and the output of the system look like this:

 and the response of the system looks like this:

which all looks normal.

Now when I am using Matlab's filter function to implement a low pass IIR filter with a cut-off frequency of 50 Hz and a Q of 1.41 then the graphs become like this:

I don't understand why the output of the system when I am using the filter function looks like this. There is a low-level curve that surrounds the 100 Hz tone and seems to follow the trend of the filter. This is causing the transfer function of the system to become messy.

I tried to implement the same low-pass filter in the frequency domain using the bilinear transform and tried the same thing. Then the results I got were much better. The response I got was what I expected: the transfer function of a second-order low pass filter.

Could someone please explain to me where I am going wrong when I am attempting the same thing with the filter function?

thank you very much Dimitris

Below is my matlab code. This is the approach that gives me the results I can't interpret. Commented out is my own implementation of the biquad while I was trying to test the result of the filter() function.

here I am attaching my matlab code:

if true
    [data, fs, nbits] = wavread(filename);

gain = 0;
f0 = 1000;
fs = 48000;
npoints = length(data);
Q = 0.707;

alpha = sin(w0)/(2*Q);
a0 =   1 + alpha;
a1 =  -2*cos(w0);
a2 =   1 - alpha;
b0 =  (1 - cos(w0)) * gainLinear / 2;
b1 =  1 - cos(w0)  * gainLinear;
b2 =  (1 - cos(w0)) * gainLinear / 2;

Ac(1)=a0/a0;
Ac(2)=a1/a0;
Ac(3)=a2/a0;
Bc(1)=b0/a0;
Bc(2)=b1/a0;
Bc(3)=b2/a0;

output = filter(Bc,Ac,input);

% x0 = 0;
% x1 = 0;
% x2 = 0;
% y1 = 0;
% y2 = 0;
% 
% output = zeros(1,length(data));
% 
% for i = 1:1:length(data)
%     x0=data(i);
%     y0 = (b0*x0 + b1*x1 + b2*x2 - a1*y1 - a2*y2)/a0;
%     y2 = y1;
%     y1 = y0;
%     x2 = x1;
%     x1 = x0;
%     
%     output(i) = y0; 
% end

NFFT = length(output); 
Y = fft(output,NFFT)*(1/fs); %multiply by period to compensate matlab scaling...
freqAxis = fs/2*linspace(0,1,NFFT/2); %construct the frequency axes
Y = 2*Y(1:NFFT/2); %keep the first half of the data, multiply by 2 to compensate for lost energy

semilogx(freqAxis,20*log10(abs(Y)),'b','lineWidth',1)
end

[buscuit]

anybody?

[Rick Rosson]

Please post your MATLAB code.

[buscuit]

hi Rick, I have appended the code in the question. Thanks