"Fred Marshall" <fmarshallx@remove_the_x.acm.org> wrote:> "Luiz Carlos" <oen_no_spam@yahoo.com.br> wrote: > > Can you point out some tool for designing FIR filters that are not > > linear phase? The ones I know only generate linear phase FIR filters! > > Well, you can use something like Octave or Matlab or Scilab to determine the > roots of the polynomial represented by the filter coefficients. [...] I don't know > of a tool that just does it in one step. There may be one.See: http://gmeteor.sourceforge.net/ which can design optimal FIR filters with arbitrary frequency responses, using a variety of optimality criteria. Cordially, Steven G. Johnson
Fast response filter?
Started by ●January 15, 2004
Reply by ●January 16, 20042004-01-16
Reply by ●January 16, 20042004-01-16
Fred Marshall wrote:> > > If you want a minimum phase FIR version of any mixed or linear phase FIR > try: > > > > -----------------------Matlab code follows--------------------------- > > function y = minphase(x) > > %minphase - minimum phase from mixed phase FIRs > > % y = minphase(x) returns the minimum phase versions of the FIRs > > % in the columns of x > > > > %Double the impulse response lengths with zeros > > x = [x; zeros(size(x))]; > > %Number of impulse responses > > m = size(x,2); > > %Size of each > > n = size(x,1); > > %The weighting function > > wn = [ones(1,m); 2*ones(n/2-1,m); ones(1,m); zeros(n/2-1,m)]; > > %The work > > y = real(ifft(exp(fft(wn.*real(ifft(log(abs(fft(x))))))))); > > %Return half the result > > y = y(1:n/2;); > > -----------------------End of code---------------------------------- > > Bob, > > Input this linear phase filter: [-1 0 2 -1 2 -1 2 0 -1] and you will find > the resulting output from above has roots that fall outside the unit circle. > So, that can't be right can it?Sure doesn't seem right.> > I found another version that looks a bit different and the results are > different. The roots fall inside the unit circle. > function y = minphase(x) > %MINPHASE minimum phase > % MINPHASE(x) returns minumum phase versions of impulse responses in x > % y = MINPHASE(x) returns the minimum phase impulse responses derived from > x. > % x = column matrix of impulse responses > m = size(x,2); > n = size(x,1); > odd = fix(rem(n,2)); > wn = [ones(1,m); 2*ones((n+odd)/2-1,m) ; ones(1-rem(n,2),m); > zeros((n+odd)/2-1,m)]; > y = real(ifft(exp(fft(wn.*real(ifft(log(abs(fft(x))))))))); > > That last expression is a mouthful isn't it!!Yeah, it's derived from the real cepstrum method that's in the Matlab rceps() function. Those lines are mine too before I found that it was necessasary to double the length first and then take half of that from the intermediate result to avoid a final result with a strange discontinuity in the middle of the IR. This is a (perhaps unwise) modification of a method I really don't understand all that well and it seems that what my modification does is sub-optimal in the sense of all the roots really being inside but is a compromise that eliminates that (perhaps necessasary) discontinuity while still keeping the result prompt.> I don't know but it seems simpler (well, easier to understand) to do > something like this: > x=...... > roots(x) > [then those roots of x that are outside the unit circle, replace with their > reciprocal real and imaginary parts] > I haven't figured out what to do with scaling.....Problem for me is that I deal with audio band FIR's much too long for roots(). Bob -- "Things should be described as simply as possible, but no simpler." A. Einstein
Reply by ●January 16, 20042004-01-16
oen_no_spam@yahoo.com.br (Luiz Carlos) writes:> Can you point out some tool for designing FIR filters that are not > linear phase? The ones I know only generate linear phase FIR filters!ScopeFIR: http://www.iowegian.com/loadfir.htm> I've made this question because sometimes we have a very small period > of time to detect one frequency, and normal (FIR, Butterworth, Cauer, > Chebichev...) band pass filters tend to have a "long" response time. > But maybe we can use another aproach, another transform that has a > faster response. > > I haven't tryed, but if we take some samples of our input data, window > these samples, padd with a lot of zeroes, and then make a DFT, we can > have fast and good precision spectrum estimator, isn't it?No. For one thing, if your signal has any noise in it, the estimate from the FFT is not optimal. Try searching for the topics "periodigram" or "Welch estimator" for more info.> I've made and interesting experience yesterday. I took a low pass > filter (Fs=8000Hz, Fpass=1000Hz, band pass ripple=0.5dB, Faten=1700Hz, > aten=45dB, Cauer, bilinear transform => fourth order filter), and fed > it whith: > -infinite to 0: signal = +1 > 0 to +infinite: signal = cos (n*2*pi*500/8000) (500Hz) > (infinite = 1000 points!!!) > Note that there is no discontinuities in this signal or it's > derivative, so it is contains only two frequencys, 0Hz until 0, and > then 500Hz.Luiz, here's a thought experiment: What would the time-domain signal look like for a spectrum consisting of only two frequencies: one at 0 Hz and one at 500 Hz? (Hint: It is NOT the signal you constructed.) If you can reach this conclusion, then you can also conclude that your signal doesn't just have 0 Hz and 500 Hz in it. Also, even if the signal did just have two frequencies in it, I'll bet you're using the PSD function or maybe just plain old "fft()" to get the spectrum, and the fft itself introduces windowing that would "splatter" your two frequencies. You gotta be careful about how you think about this stuff.> But the filter "created" a "spectral noise" to transit from 0Hz to > 500Hz, even though the two frequencys are in it passband!No, the filter creates no such thing, not in theory at least. Filters are linear systems, one of the properties of which is that they do not introduce new frequencies. Now of course if you screw up in your *implementation* of the filter, like not scaling it properly so that you get overflows, you will get new frequencies. But that's not the filter's fault - that's yours!> And this "deformation" is proportional to the response time. > So I'm asking myself what those filters in the DA converter are doing > with my "pure" audio???Filtering.> Sometime ago, someoneMe.> asked why we are using so high sampling > frequencys in the new digital audio standards. Maybe here is the > response. Higher sampling frequencys implies in more transition band > for the low pass filters (DA stage), with less taps and then, less > distortion.No, but nice try Luiz. It shows you are thinking about things, and that is a good thing. :) Keep the thoughts and questions coming. -- % Randy Yates % "Bird, on the wing, %% Fuquay-Varina, NC % goes floating by %%% 919-577-9882 % but there's a teardrop in his eye..." %%%% <yates@ieee.org> % 'One Summer Dream', *Face The Music*, ELO
Reply by ●January 19, 20042004-01-19
Thanks to everybody!
Step by step, the easy ones first:
1- Randy, how can I use ScopeFIR to generate a Minimum Phase Filter? I
could not find!
2- About the transition band:
Maximum frequency to be maintained = 20kHz.
Fs=48kHz, that means our transition band is 4kHz wide to avoid
aliasing.
Fs=96kHz, now we have a 28kHz wide transiton band. Right?
3- About my signal (Jerry, Robert and Randy):
Well, I did it's FFT.
A 16 point FFT, input = {1, 1, 1, 1, 1, 1, 1, 1, 1, 0, -1, 0, 1, 0,
-1, 0}.
(1 followed by Fs/4).
Modulus result= {8, 5.4, 0, 3.9, 4, 1.7, 2.1, 0}.
So, really there are frequencies higher than Fs/4!!!! Very, very
surprising to me.
I can feel why a transition from one value to another (a
discontinuity), or a transition from zero to a sine (a knee, first
derivative discontinuity), both are spectrally spread. But my
example... I simply just can't see where are these high frequency
components coming from! I always relate the spectral distribuition of
a signal to the speed it changes, ..., looks like there is something I
am missing.
Any help will be very appreciated.
Luiz Carlos.
Reply by ●January 19, 20042004-01-19
Luiz Carlos wrote: ...> I can feel why a transition from one value to another (a > discontinuity), or a transition from zero to a sine (a knee, first > derivative discontinuity), both are spectrally spread.OK think about superposing the components of your signal at negative and positive time. Each has the discontinuity you understand and the spreading spectrum caused by that discontinuity. There is no reason to assume that the transients of the one cancel the transients of the other. Experiment shows that they do not.> But my > example... I simply just can't see where are these high frequency > components coming from!The "Argument from Disbelief" is dangerously seductive -- creationists use it often -- but not logically sound.> I always relate the spectral distribuition of > a signal to the speed it changes, ..., looks like there is something I > am missing.Attack the problem from the other end. What combination of sines and cosines of various frequencies, existing throughout time, would you need to add together to create your test waveform? Those are its components; anything fewer won't do. (You can find them by FFTing the waveform you created.)> Any help will be very appreciated.I hope you have something to appreciate. If not, someone will try again.> Luiz Carlos.Jerry -- Engineering is the art of making what you want from things you can get. �����������������������������������������������������������������������
Reply by ●January 19, 20042004-01-19
Luiz Carlos wrote:> I can feel why a transition from one value to another (a > discontinuity), or a transition from zero to a sine (a knee, > first derivative discontinuity), both are spectrally spread. But > my example......has a discontinuous second derivative. Martin
Reply by ●January 20, 20042004-01-20
>...has a discontinuous second derivative.Yes Martin, and as a sum of sinusoids, all derivatives must be continuous, isn't it? Please read the rest. Jerry, my confusion comes from a very long time ago! When I learned about the sampling theorem, the Nyquisq (Shannon) frequency was explained intuitively using the carriage wheel presented in a movie. If it spins too fast you see it spinning backwards (aliasing). In other words, you must sample at least two times faster than the faster "thing" you are capturing. So, when I look at a signal, I do my "visual FFT", observing how it varies on time, not in frequency, because frequency are consequence of variation in time. Discontinuities, knees, and very high first derivatives are easy seen. But as Martin pointed out, the second derivative must be continuous too. And looking more carefully, all derivatives must be continuous. Now I'm disappointed, because my mental FFT algorithm has a flaw! :)> I hope you have something to appreciate. If not, someone will try again.I can't say I'm happy to kwow I'm wrong. But I really prefer to know that I'm wrong! It's good to learn somethig new, even though it was all the time in front of my eyes. Jerry, I did what you suggested. The FFT, etc. But I needed to understand what was happening. Thankyou, Luiz Carlos
Reply by ●January 20, 20042004-01-20
"Bob Cain" <arcane@arcanemethods.com> a �crit dans le message de news:40088F3E.5CE04EED@arcanemethods.com...> > > Fred Marshall wrote: > > > > > If you want a minimum phase FIR version of any mixed or linear phaseFIR> > try: > > > > > > -----------------------Matlab code follows--------------------------- > > > function y = minphase(x) > > > %minphase - minimum phase from mixed phase FIRs > > > % y = minphase(x) returns the minimum phase versions of the FIRs > > > % in the columns of x > > > > > > %Double the impulse response lengths with zeros > > > x = [x; zeros(size(x))]; > > > %Number of impulse responses > > > m = size(x,2); > > > %Size of each > > > n = size(x,1); > > > %The weighting function > > > wn = [ones(1,m); 2*ones(n/2-1,m); ones(1,m); zeros(n/2-1,m)]; > > > %The work > > > y = real(ifft(exp(fft(wn.*real(ifft(log(abs(fft(x))))))))); > > > %Return half the result > > > y = y(1:n/2;); > > > -----------------------End of code---------------------------------- > > > > Bob, > > > > Input this linear phase filter: [-1 0 2 -1 2 -1 2 0 -1] and you willfind> > the resulting output from above has roots that fall outside the unitcircle.> > So, that can't be right can it? > > Sure doesn't seem right.why ? A linear phase filter doesn't imply that it's minimum phase. I think they are always non minimum phase but I don't remember why. Patrick
Reply by ●January 20, 20042004-01-20
Patrick wrote:> > > > Bob, > > > > > > Input this linear phase filter: [-1 0 2 -1 2 -1 2 0 -1] and you will > find > > > the resulting output from above has roots that fall outside the unit > circle. > > > So, that can't be right can it? > > > > Sure doesn't seem right. > > why ? A linear phase filter doesn't imply that it's minimum phase. I think > they are always non minimum phase but I don't remember why. >The intent of the calculation is to produce a minimum phase filter from a mixed or linear phase one. My modification, while seeming to give a better time domain function, misses the mark of being truly minimum phase. Close though. :-) Bob -- "Things should be described as simply as possible, but no simpler." A. Einstein
Reply by ●January 20, 20042004-01-20
"Patrick" <busymanfr@yahoo.fr> wrote in message news:bujei9$l2i@aix4.segi.ulg.ac.be...> > "Bob Cain" <arcane@arcanemethods.com> a �crit dans le message de > news:40088F3E.5CE04EED@arcanemethods.com... > > > > > > Fred Marshall wrote: > > > > > > > If you want a minimum phase FIR version of any mixed or linear phase > FIR > > > try: > > > > > > > > -----------------------Matlab codefollows---------------------------> > > > function y = minphase(x) > > > > %minphase - minimum phase from mixed phase FIRs > > > > % y = minphase(x) returns the minimum phase versions of the FIRs > > > > % in the columns of x > > > > > > > > %Double the impulse response lengths with zeros > > > > x = [x; zeros(size(x))]; > > > > %Number of impulse responses > > > > m = size(x,2); > > > > %Size of each > > > > n = size(x,1); > > > > %The weighting function > > > > wn = [ones(1,m); 2*ones(n/2-1,m); ones(1,m); zeros(n/2-1,m)]; > > > > %The work > > > > y = real(ifft(exp(fft(wn.*real(ifft(log(abs(fft(x))))))))); > > > > %Return half the result > > > > y = y(1:n/2;); > > > > -----------------------End of code---------------------------------- > > > > > > Bob, > > > > > > Input this linear phase filter: [-1 0 2 -1 2 -1 2 0 -1] and you will > find > > > the resulting output from above has roots that fall outside the unit > circle. > > > So, that can't be right can it? > > > > Sure doesn't seem right. > > why ? A linear phase filter doesn't imply that it's minimum phase. Ithink> they are always non minimum phase but I don't remember why. >Patrick, ??? Being linear phase means that it's definitely not minimum phase. I don't think either Bob or I suggested anything different. - The reason why is because a linear phase filter has reciprocal pairs of zeros inside and outside the unit circle. - If the linear phase filter has zeros that are complex, then they are paired with their complex conjugates - making each complex zero instance a quad and each real zero instance a pair. - For minimum phase, all the zeros have to be inside the unit circle. So the two characteristics are mutually exclusive. - If any complex conjugate zero pair is replaced with its reciprocal, the magnitude response is unchanged. - If all the zeros are outside the unit circle then it's maximum phase. - So, a linear phase filter is a combination of minimum and maximum phase. - And, there are a multiplicity of filters (of the same order) with the same magnitude response that are had by moving zeros inside or outside the unit circle. - Note that a linear phase filter is really a cascade of a minimum phase and a maximum phase so that its magnitude response is like the square of the two and not the same as either one. Another way to look at this is that you can get a minimum phase filter with "similar" magnitude response by taking the minimum phase square root of a linear phase filter if it's designed to have all positive mangitude response to begin with (Schuessler). Fred
