Mark Borgerding <mark@borgerding.net> wrote in message news:<fCLkd.143680$5v2.39391@fe2.columbus.rr.com>...
> Rune Allnor wrote:
> > Mark Borgerding <mark@borgerding.net> wrote in message news:<9gqkd.1761$DI5.1733@fe1.columbus.rr.com>...
> > 
> > If you read the exchange of opinion between Randy and me, you will find 
> > that one opinion goes as "Non-FIR filters filters can not be implemented 
> > in frequency domain" while the other goes as "yes they can, and here's 
> > how it can be done". Regardless of practical issues, I think the two 
> > points of views are mutually exclusive, and I find it hard to see that 
> > both can be right.
> 
> Randy says one can only use fast convolution for FIR filters.

Well... it's a matter of interpreting the original question. I agree 
if you say "one can only use fast convolution for finite amounts 
of data". In that sense, any practical data set is finite, and an IIR 
filter does not really yield an infinite impulse response and can be 
said to be FIR. My interpretation of the question is whether recursive 
filters (I'm rephrasing here, to avoid splitting hairs over the 
word "infinite") can be implemented in frequency domain.

True, the price paid for using the frequency domain method for IIRs 
is the risk of wrap-around effects, but with care you can implement 
the IIR as a frequency domain solution. Whether that's an easy or 
smart thing to do, is a completely different matter. 

> You show how you can implement IIR filtering with fast convolution ... 
> by turning that IIR filter into an FIR filter.

Yep. That's a common way of doing things. The main benefit when 
compared to the "IIR way" is that you have complete control of end 
transients. You do need to compute the impulse response, though, from 
the IIR filter coefficients. How would you du that? Some use the 
unit impulse sequence and implement the difference equation.
That way one can decide the "effective N" on the fly, by monitoring 
the magnitude of the response samples. It would require a "sufficiently 
large" data buffer being available in the first place, to store the 
coefficients. If that's not available, you would have to estimate 
the impulse response twice: The first time to estimate the number of 
"effective" coefficients in the impulse response; the second time 
to compute the coefficients before storing them after having allocated
the memory space. 

An alternative would be to try

h=real(ifft(fft(b,N)./fft(a,N)));

where a and b are the filter coefficient vectors and N is "sufficiently 
large" to contain the significant parts of the impulse response. 
Again, the problem is that N must be specified in advance. I don't say 
that this is a particularly smart way of doing things, but it's 
_possible_ to do it like this.

> So who's right? Both of you.
> Although Rune gets bonus points for creative solutions.

Thanks for the credit, but I picked it up somewhere on comp.dsp. 
I can't remember who posted the trick, though.
 
> Myself, I prefer to do the IIR to FIR conversion explicitly in the time 
> domain, rather than guessing at transient lengths

Agreed. That's the way one would do it in practice. Again, my 
interpretation of the question was whether pure frequency domain 
operations are _possible_.

> The question we should be asking is, "why use IIR filters if you can use 
> fast convolution?"

There was a thread here, some time ago, where the poster tried to 
do some processing in frequency domain and convert back to time 
domain and play the result on the speaker. The resulting sound 
wasn't very good, and I think glitches at the joints between blocks 
was at least part of the problem. Frequency domain solutions might 
work in off-line applications. In on-line applications you might get 
into other kinds of trouble that time-domain implementations avoid.
 
> The OP may have a legitimate reason, but it is good to keep sight of FIR 
> design methods in the process.  A (The?) reason for using IIR filters is 
> to get better filter shape than a direct-implemented FIR of equal cost 
> could provide.

There are other reasons. A student of mine did some data processing 
where it turned out that control of passband ripple was essential for
getting good results. The trick only worked when he used Butterworth 
or Cheb 2 filters. As far as I know, he used the IIR implementation. 
I don't know if things would have been messed up if he had used a FIR 
approximation (actually, I doubt that a FIR approximation would have 
made much of a difference).

Rune

Mark Borgerding <mark@borgerding.net> writes:
> [...]
> Randy says one can only use fast convolution for FIR filters.
> 
> You show how you can implement IIR filtering with fast convolution
> ... by turning that IIR filter into an FIR filter.

You're attempting to show some grace, Mark - a personal quality of
the highest order.

The hard fact is, "finite != infinite," just as much as "is != (is
not)." Thus the finite DFT (as defined by Oppenheim et al. in "Signals
and Systems," and which is by definition the tool used in
frequency-domain filtering as defined in P&M) can NEVER be used to
implement an infinite impulse response.

And to put a stop to the "practical implementation" argument, I was
speaking analytically, not practically.

(Please side-step the emotional power in this message, Mark et al.,
and allow it to reach its intended target.)
-- 
Randy Yates
Sony Ericsson Mobile Communications
Research Triangle Park, NC, USA
randy.yates@sonyericsson.com, 919-472-1124

Rune Allnor wrote:
> Mark Borgerding <mark@borgerding.net> wrote in message news:<9gqkd.1761$DI5.1733@fe1.columbus.rr.com>...
> 
> If you read the exchange of opinion between Randy and me, you will find 
> that one opinion goes as "Non-FIR filters filters can not be implemented 
> in frequency domain" while the other goes as "yes they can, and here's 
> how it can be done". Regardless of practical issues, I think the two 
> points of views are mutually exclusive, and I find it hard to see that 
> both can be right.

Randy says one can only use fast convolution for FIR filters.

You show how you can implement IIR filtering with fast convolution ... 
by turning that IIR filter into an FIR filter.

So who's right? Both of you.
Although Rune gets bonus points for creative solutions.


Myself, I prefer to do the IIR to FIR conversion explicitly in the time 
domain, rather than guessing at transient lengths


The question we should be asking is, "why use IIR filters if you can use 
fast convolution?"

The OP may have a legitimate reason, but it is good to keep sight of FIR 
design methods in the process.  A (The?) reason for using IIR filters is 
to get better filter shape than a direct-implemented FIR of equal cost 
could provide.


Example: Butterworth(5) vs. Remez(32)

A 5th order butterworth (IIR) filter can be expressed (with error < 
-90dB) by an impulse response of length 34.

A Remez filter order 32 (33 coeffs), has a shorter impulse response, a 
better frequency response, and linear phase.

But don't take *my* word for it ...

% Here's a little Octave/Matlab code to demonstrate IIR to FIR
% truncation and compare it to similar length Remez filter
% Copyright Mark Borgerding, 2004
% Licensed under a BSD License
% see http://www.fsf.org/licenses/info/BSD_3Clause.html

[bbi,abi]=butter(5,.5);
bbf = filter(bbi,abi,[1 zeros(1,100)]);
bbfpow = sum(bbf.^2);
trunc_pow_loss = 10*log10( 1 - cumsum(bbf.^2)/bbfpow);
coefs_needed = min( find( trunc_pow_loss < -90) )
% coefs_needed is 34

bremez = remez(32,[0 .3 .8 1],[1 1 0 0]);

figure(1);
title('original IIR Butterworth');
freqz(bbi,abi);

figure(2);
title('FIR approx of Butterworth');
freqz(bbf(1:coefs_needed) );

figure(3);
title('Remez filter of similar cost to FIR approx');
freqz(bremez);


Cheers,
   Mark


<SHAMELESS_PLUG>
If you like what you see above, and would like to inquire how I could 
help your project, contact markb@3dB-Labs.com for more info.

If you don't like what you see above, and think I can't even spell DSP, 
then I claim this post was written a roaming gang of rabid monkeys.
Dang monkeys.
</SHAMELESS_PLUG>

Mark Borgerding <mark@borgerding.net> wrote in message news:<9gqkd.1761$DI5.1733@fe1.columbus.rr.com>...
> Rune Allnor wrote:
> > Randy Yates <yates@ieee.org> wrote in message news:<actt9icg.fsf@ieee.org>...
>  ...
> >>Note that frequency-domain filtering is only useful for FIR filters, i.e.,
> >>filters in which the a vector above is 1. Butterworth filters are not
> >>FIR filters, so you won't be able to use frequency-domain filtering on them.
> > 
> > 
> > Yes you are. With the filter coefficients grouped as a and b above, 
> > the frequency domain IIR filter is implemented as
> > 
> > X=fft(x,N);           % DFT of input signal
> > A=fft(a,N);
> > B=fft(b,N); 
> > Y= X.*B./A;
> > y= real(ifft(y));
> > 
> > where the FFT length N is large enough to contain both the input signal 
> > and startup and end transients. The only difficult part here is to 
> > find a suitable N. Which can be hard enough, though.
> > 
> > Rune
> 
> Rune,
> This is really just an approximation of an IIR filter by taking the 
> first (most significant) portion of its impulse response and truncating 
> it, thus making it a finite impulse response.

Uh, yes? 

Of course, any implementation that works on a finite number of samples 
can only be an approximation of a true IIR. I don't know of many 
practical implementations that are truly infinite.

> With a stable IIR filter, this approximation can be made close enough 
> for any practical purposes.

Which, of course, is the interesting bit.
 
> I'd say Rune and Randy are both right, depending on the viewpoint.

If you read the exchange of opinion between Randy and me, you will find 
that one opinion goes as "Non-FIR filters filters can not be implemented 
in frequency domain" while the other goes as "yes they can, and here's 
how it can be done". Regardless of practical issues, I think the two 
points of views are mutually exclusive, and I find it hard to see that 
both can be right.

> In the strict sense, it is imposible to implement an infinite impulse 
> response with circular convolution.  

I'd go a bit further and say "it is imposible to implement an infinite 
impulse response" regardless of whether one uses circular convolution 
or not. Of course,  in the case of the FIR one has a finite end 
transient, so the practical aspects of determining the N above become 
easier (N = length(x) + 2*length(b)-2;). In the IIR case, one needs to 
find a "effective length" of the impulse response, and decide at what 
time the impulse response no longer makes significant contributions to 
the start/end transient. This needs not be an easy task, but it is 
possible to do within the constraints of practical limitations to IIR 
systems.

Now, the question was whether it is _possible_ to implement a Butterworth 
filter (I'll take the liberty to generalize this to a 'recursive filter') 
in frequency domain. If you want to go hypothetical and discuss the 
infinite time case, use the 'proper' DFT (with infinite summation ranges) 
instead of the FFT to compute the spectra A and B above. The theory works 
perfectly well in frequency domain. You only need the filter 
coefficients, no time-domain data at all.   

> In practice, you can get as close 
> as you want (excepting instability/marginal stability).

Agreed.

Rune

Rune Allnor wrote:
> Randy Yates <yates@ieee.org> wrote in message news:<actt9icg.fsf@ieee.org>...
...
>>Note that frequency-domain filtering is only useful for FIR filters, i.e.,
>>filters in which the a vector above is 1. Butterworth filters are not
>>FIR filters, so you won't be able to use frequency-domain filtering on them.
> 
> 
> Yes you are. With the filter coefficients grouped as a and b above, 
> the frequency domain IIR filter is implemented as
> 
> X=fft(x,N);           % DFT of input signal
> A=fft(a,N);
> B=fft(b,N); 
> Y= X.*B./A;
> y= real(ifft(y));
> 
> where the FFT length N is large enough to contain both the input signal 
> and startup and end transients. The only difficult part here is to 
> find a suitable N. Which can be hard enough, though.
> 
> Rune

Rune,
This is really just an approximation of an IIR filter by taking the 
first (most significant) portion of its impulse response and truncating 
it, thus making it a finite impulse response.

With a stable IIR filter, this approximation can be made close enough 
for any practical purposes.


I'd say Rune and Randy are both right, depending on the viewpoint.
In the strict sense, it is imposible to implement an infinite impulse 
response with circular convolution.  In practice, you can get as close 
as you want (excepting instability/marginal stability).


-- Mark

snsundar@olemiss.edu (sundar) writes:

> yes even i do , thank you so much. 

Sure thing, sundar.

> Randy, i have one more question, after i do a 
> hilbert(y)
> ang=angle(y);
> plot(ang(:,1))
> 
> this plot that iam getting has a pattern for a signal and is not
> random for noise. Can you tell me why, am i not getting a random phase
> angle distribution for that one sec data for noise. I read some of the
> old posts and think that i might have something to do with my
> filtering, and since my filtering widht is very narrow, all my points
> are sort of correlated and hence lose their random behavior and show a
> pattern.

I think I agree with the "old posters." Here's one way to look at it:

The output power-spectral density Sy(w) of a signal x passed throw a filter 
H(w) is 

  Sy(w) = Sx(w) * |H(w)|^2,

where Sx(w) denotes the input signal's PSD. Also, the Wiener-Khinchine
theorem states that the power spectrum of a signal Sy(w) is the Fourier
transform of its autocorrelation function Ry(k), or said from the
other direction, the autocorrelation Ry(k) of a signal is the *inverse*
Fourier transform of its PSD,

  Ry(k) = T_s/(2*pi) \int_{\pi/T_s}^{+\pi/T_s} Sy(w) e^{+j*w*k*T_s) dw.

So if you have an H(w) that is very narrow, the resulting Sy(w) is
also very narrow, and the resulting inverse FT will be "wide," i.e.,
non-zero at a lot of points, indicating there is a lot of correlation
in the signal. 

A much easier way to look at it is just as you've begun: a narrow
filter produces a sine way.

> i tried 
> angun=unwrap(ang);   % unwrap my phase
> phi=angun(:,1)-2*pi*freq*(0:1/fs:1)'; % remove my carrier freq
> plot(phi)
> phi1=atan2(sin(phi1),cos(phi1));  % convert from -pi to pi
> plot(phi1)
> 
> This again leads me nowhere.

Yeah, I don't doubt it.

> Using the hilbert transform how can i differentiate between signal and
> noise.

I don't know. Who said it could be used that way? 

> thanks so much once again...u could mail me and post the reply since
> that way it wud  be faster for us to communicate.
> 
> Thanks so much
> 
> Regards
> 
> Sundar    snsundar@olemiss.edu

My Dad's family is from Taylorsville, MS - a small berg between
Laurel and Jackson. Several of my uncles and cousins have attended
Ole Miss. I'm jealous! 

Presuming you're Indian, which part are you from? I'm going to
Bangalore the end of the month to visit my fiancee. 
-- 
Randy Yates
Sony Ericsson Mobile Communications
Research Triangle Park, NC, USA
randy.yates@sonyericsson.com, 919-472-1124

yes even i do , thank you so much. 

Randy, i have one more question, after i do a 
hilbert(y)
ang=angle(y);
plot(ang(:,1))

this plot that iam getting has a pattern for a signal and is not
random for noise. Can you tell me why, am i not getting a random phase
angle distribution for that one sec data for noise. I read some of the
old posts and think that i might have something to do with my
filtering, and since my filtering widht is very narrow, all my points
are sort of correlated and hence lose their random behavior and show a
pattern.

i tried 
angun=unwrap(ang);   % unwrap my phase
phi=angun(:,1)-2*pi*freq*(0:1/fs:1)'; % remove my carrier freq
plot(phi)
phi1=atan2(sin(phi1),cos(phi1));  % convert from -pi to pi
plot(phi1)

This again leads me nowhere.

Using the hilbert transform how can i differentiate between signal and
noise.

thanks so much once again...u could mail me and post the reply since
that way it wud  be faster for us to communicate.

Thanks so much

Regards

Sundar    snsundar@olemiss.edu

snsundar@olemiss.edu (sundar) writes:
> [...]
> ffty=fft(y); % at this step here i should be able to c my peak at
> 380hz but dont..

I do. Are you typing, e.g., 

plot(abs(ffty(:,1)));

? (you have 60 columns so you could use 1 to 60 for the second index).
-- 
%  Randy Yates                  % "Watching all the days go by...    
%% Fuquay-Varina, NC            %  Who are you and who am I?"
%%% 919-577-9882                % 'Mission (A World Record)', 
%%%% <yates@ieee.org>           % *A New World Record*, ELO
http://home.earthlink.net/~yatescr

Hi Randy and all others, 

thank you for your help and timely replies. I am really desperate to
get this done asap, my graduation depends on this.

yes,i understand this, but my binwidht is one. 
> Well, it depends on what you mean by "one peak". The windowing characteristic
> of the FFT will smear a peak, but the maximum should still be at 380 Hz. Note that depending on the sample rate and FFT length, 380 Hz may or may not correspond
to a specific FFT bin, so you may have more than one bin at the peak
value.
or right now, while using Matlab, the way you're going (using
filter())

> sundar, can you post your Matlab code and your data to a web page somewhere?
That way we can see exactly what's going on.
Well i dont have a website..let me trying copying the code here...

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% initial parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% sampling rate for 60 sec
fs=2100

%fft_size
fft_size=2100

%binwidht
bwidht=fs/fft_size

%desired frequency
freq=380

%size of filter to apply
df=1

% apply a butterworth filter across the desired freqs.

w1 = (freq-df/2)/fs
w2 = (freq+df/2)/fs

% this gives me w1 at 379.5 & w2 at 380.5

% order of the butterworth filter

n=2
Wn=[w1 w2]

[b,a]=butter(n,Wn);

y=filter(b,a,data); % data(2100,60)

ffty=fft(y); % at this step here i should be able to c my peak at
380hz but dont..

% hilbert tranform remove f<0, 2*f>0;

for i=1:1:fft_size
	if i<=fft_size/2
	nfftdata(i,:)=ffty(i,:).*2;
else
	nfftdata(i,:)=0;
end
end

hildata=ifft(nfftdata,fft_size);

ang=angle(hildata);

% hence this will give me 2100 angle points for my first second data.

I need to prove that my angle should be random for noise while it
should not be random for signals. The second part works fine with a
hann windowing but i get a pattern for noise which i dont think is
logical. while with this filter function i get the same angle plot for
signal and noise.

Thanks once again.

Regards

Sundar

snsundar@olemiss.edu (sundar) writes:

> Randy Yates <yates@ieee.org> wrote in message news:<actt9icg.fsf@ieee.org>...
> > snsundar@olemiss.edu (sundar) writes:
> Dear Randy 
> 
> Thanks once again for your advice. I tried using the filter function,
> but iam not sure if i got what is right. Now, ok i have a signal at
> 380 Hz, i use a butterworth band pass filter from 379.5 to 380.5. Get
> [b.a], then use it on my freq domain data. Now shouldnt my fft plot
> show only one peak at 380 hz.

Well, it depends on what you mean by "one peak". The windowing characteristic
of the FFT will smear a peak, but the maximum should still be at 380 Hz. Note
that depending on the sample rate and FFT length, 380 Hz may or may not correspond
to a specific FFT bin, so you may have more than one bin at the peak value. 

> I think iam getting lost in terminology (freq domain filtering etc. )
> above mentioned is my objective what do you think is the best method
> of doing it, 

For right now, while using Matlab, the way you're going (using filter())
is just fine. 

It sounds like what you need to do filter a signal (whether you use a 
time-domain filter or frequency-domain filter is an implementation issue
which you probably don't need to consider right now) and look at the 
resulting response in the frequency domain. So I think you're just saying
you want to use the FFT to view the spectrum of your filter's output, which
is a valid thing to do.

> i tried rectangular window and hanning, they dont seem to
> help.

sundar, can you post your Matlab code and your data to a web page somewhere?
That way we can see exactly what's going on.
-- 
Randy Yates
Sony Ericsson Mobile Communications
Research Triangle Park, NC, USA
randy.yates@sonyericsson.com, 919-472-1124