I'm not a Spice user, but all you need to do is take a normal FFT and apply a tilt to it, then do an IFFT (if you want to see the time waveform). Doesn't Spice allow you to do this?
The tilt should be +6 dB/octave for a normal derivative (see my prior post) or in your case some fraction of that. Note that the raw FFT is linear in frequency and the "tilt" I am talking about is linear in *log* frequency, so you need a little exponential math here to generate the proper shape.
Best regards,
Bob Masta DAQARTA v6.00 Data AcQuisition And Real-Time Analysis
formatting link
Scope, Spectrum, Spectrogram, Sound Level Meter Frequency Counter, FREE Signal Generator Pitch Track, Pitch-to-MIDI Science with your sound card!
The usefulness of windowing isn't restricted to the noiseless CW case. Unless your signal is actually periodic with the same period as the DFT, failure to window will cause every bin to smear out all across the interval. You may not notice it with white noise, because the spectrum stays flat in that case, but there are other kinds of noise where you certainly will.
If your signal is a pulse, so that it starts from zero somewhere in the interval and decays back to zero by the end, then you effectively are in the periodic case already. If not, windowing is good medicine.
Cheers
Phil Hobbs
--
Dr Philip C D Hobbs
Principal
ElectroOptical Innovations
55 Orchard Rd
Briarcliff Manor NY 10510
845-480-2058
email: hobbs (atsign) electrooptical (period) net
http://electrooptical.net
The best you can do for differentiation in a DFT is to take the central finite difference, but the DFT of that operation is not a ramp, linear, logarithmic, or exponential.
With an appropriately windowed transform, you can make a DFT look like samples of the continuous-time Fourier transform of a band-limited continuous time function.
Cheers
Phil Hobbs
--
Dr Philip C D Hobbs
Principal
ElectroOptical Innovations
55 Orchard Rd
Briarcliff Manor NY 10510
845-480-2058
email: hobbs (atsign) electrooptical (period) net
http://electrooptical.net
I don't know what you mean by "best" re: the central finite difference in view of the need for using finite discrete transforms and/or windowed data. As far as I know, a minimax FIR filter designed differentiator approximates an imaginary ramp in frequency *up to* some cutoff frequency - which is necessary to avoid aliasing I do believe. Another approach to getting this is a Taylor approximation to the central finite difference of the first derivative which apparently comes in closed form (Khan and Ohba, 1999). I believe these approaches are common practice in a case like this.
Here is an example of a length 51 differentiator of the minimax type done using the Parks-McClellan program: [2.83E-05 -7.68E-05 1.34E-04 -2.24E-04 3.53E-04 -5.32E-04 7.73E-04 -1.09E-03 1.50E-03 -2.01E-03 2.66E-03 -3.46E-03 4.44E-03
Thanks, that's pretty. I was perhaps a bit unclear in saying that the first finite difference is the best you can do, but the point at issue was what happens at frequencies near +-f_sample/2. The first finite difference transforms to a sine wave with a 1/2 sample delay.
Optimal differentiators can be looked at as a variant of windowing a ramp, by using linear combinations of the first M finite differences to make an M+1 order approximation to the derivative over some frequency band, but they always roll off at high frequency, which is a key point.
Bob M. seems to be suggesting using a real ramp with no high-frequency cutoff, which will give what the compiler manuals describe as "unexpected results".
Cheers
Phil Hobbs
--
Dr Philip C D Hobbs
Principal
ElectroOptical Innovations
55 Orchard Rd
Briarcliff Manor NY 10510
845-480-2058
email: hobbs (atsign) electrooptical (period) net
http://electrooptical.net
Yes. I'm glad you didn't say cosine. It's an imaginary sine wave at the lowest frequency. So, it has one period. As such, it goes from positive to negative at fs/2. And, it's "ramp like" near f=0 as sin(x) = x for x < whatever small number (0.25 is 1% error) So, it's a reasonable differentiator for frequencies below whatever limit you like there. Not so good at pi/2 or fs/4!!
The "imaginary" part seems to be getting lost in the "ramp" terminology. Note that the FIR filter I provided is real and purely odd in time - thus odd and purely imaginary in frequency (if centered at t=0). So, any of these is an approximation to "jw" - an imaginary ramp.
You will still get a better behaved result (ie one which is mainly due to the core transient with more predictable and well controlled artefacts if you put that in the middle of the window and tile it with a phantom unmeasured vertical axis reflection time shifted (or rather do a transform that uses that implied symmetry).
Otherwise you will have the result of your measured transient combined with the very sharp step down at the end of the range and associated high frequency components that may not be in your measured data.
Windowing to avoid sharp edge discontinuity artefacts is usually worthwhile - and for sampled data that might not be truly bandlimited and subject to aliassing it can make sense to preconvolve the raw data with a gridding function to get something that behaves much more like the ideal DFT. Prolate spheroidal Bessel functions are amongst the simplest which give real benefit in this application see for example:
formatting link
This method requires a multiplicative correction after the FFT but gives a final result much closer to the DFT.
What you say here about power of two transforms is no longer true. FFTW and other fast low prime transform kernels are now faster than naieve power of two implementations for some favourable sizes and competitive for many more. In some cases due to associative cache problems on certain processors exact powers of two are slower than they should be!
The pattern of FFT radix for which non power of two beats zero padding to the next power of two for my local FFT implementation is:
Taking only the ones where there is a clear win. Quite a lot of others are broadly competitive with the longer highly optimised 2^N transform.
3072 and 3125 are nearly 2x faster than an optimised radix-8 4096 and 3x faster than the naive radix-2 implementation for instance. An optimised split radix-16 transform would shave roughly another 10% off the 4096 time but would not change the outcome.
I wonder if there isn't some misunderstanding here? The notion of zero padding to achieve a power of two seems an archaic concern - as you well point out. However, zero padding to accommodate circular convolution is common - much more common I should think.
Nonetheless, you've pointed out some things that are new info for me. Thanks.
Interesting. To clarify, you seem to be using 'DFT' to mean a sampled-time Fourier transform of infinite length, is that right? IME that's usually used to denote the circular version for which the FFT is a fast algorithm.
I'm also less sanguine about resampling. In the 1-D case at least, the orthogonal functions for the actual sampling grid used may be quite ill-conditioned, and resampling without taking account of any special features of the eigenfunctions can lead to tears.
Cheers
Phil Hobbs
--
Dr Philip C D Hobbs
Principal
ElectroOptical Innovations
55 Orchard Rd
Briarcliff Manor NY 10510
845-480-2058
email: hobbs (atsign) electrooptical (period) net
http://electrooptical.net
Zero padding a transient without windowing does zilch to reduce spectral leakage--all it does is to evaluate the exact same transform at different sample frequencies. You can see that from the definition of the DFT (circular).
Cheers
Phil Hobbs
--
Dr Philip C D Hobbs
Principal
ElectroOptical Innovations
55 Orchard Rd
Briarcliff Manor NY 10510
845-480-2058
email: hobbs (atsign) electrooptical (period) net
http://electrooptical.net
I should add that leaving space to avoid wraparound in the circular convolution (as Fred mentioned) doesn't suffer this problem.
Cheers
Phil Hobbs
--
Dr Philip C D Hobbs
Principal
ElectroOptical Innovations
55 Orchard Rd
Briarcliff Manor NY 10510
845-480-2058
email: hobbs (atsign) electrooptical (period) net
http://electrooptical.net
Yes. I am using DFT here in the context of the idealised FT of infinite length that is zero outside the region of interest. The one that you would get by summing the discrete Fourier equations without the implicit FFT assumption of tiled or mirror boundary periodicity.
This becomes important in classical radio astronomy aperture synthesis when there are bright point sources in the field of view with grating ring artefacts from the telescope baseline configuration that cross the reconstructed image boundary!
Sorry for any confusion.
The neat thing about the prolate spheroidal Bessel functions is that they are the Eigenfunctions of bandlimiting an FT to some finite width and then taking its FT again. This means that when you regrid data using one of them all components are accurately represented. The adhoc truncated Gaussian, Kaiser-Bessel or truncated Gaussian*Sinc functions commonly used have weaker alias rejection and/or in band distortions that affect high dynamic range data.
I don't agree with everything in this paper, but it provides an introduction and comparison of several popular FT gridding functions. The older papers are not online.
formatting link
There is another even more sophisticated family of gridding functions that sacrifice a guard band around the edge of the transform to obtain an extremely high fidelity FT of interferometer and MRI data. Basically you can't trust the edges or the corners because some leakage always occurs and these methods formalise the error bound in the trusted zone.
Remember, this is in the context of converting a step response into an impulse response. The step response waveform of a perfect system is a flat line at unity... there is no step-down.
In a real system that has an overall low-pass response, the step response starts below unity and rises up to unity (maybe with some overshoot and damped oscillation) and then stays at unity until the end of the FFT period. We don't want to taper that down to zero... staying at unity is the desired result.
It's true that there is a conceptual step-down that is off-screen, to prepare for the next rising edge. But that has no effect on the desired step response.
In the case of overall high-pass systems, the step response droops down from the initial peak, and will be somewhere less than unity by the end of the frame. Although this somehow seems more worrisome than the low-pass case, in fact it is not really a problem for this approach. The fact that the step response is not satisfactorily completed means that the low end of the frequency response is not fully captured. That's all stuff that is between the 0th (DC) and 1st FFT lines... a longer FFT is needed to resolve that.
Thanks for the FFTW reference. I am gratified to see that there are finally decent non-proprietary FFT routines available for the public to use. When I got started in this adventure, back in the days of the IBM PC and XT, it was "common knowledge" that you couldn't do real-time spectral analysis on those systems without a dedicated custom add-in card, but that "someday soon" we could hope to do it all on the main CPU. I realized that I had been hearing this line for years already, in reference to the Apple, Commodore, and Sinclair (etc). It dawned on me that maybe nobody had really tried to optimize an FFT for this purpose, since all the references I could find were obviously not well optimized. They were written in FORTRAN or BASIC, computed trig functions on the fly in every iteration, and did a few other obviously redundant things. Coding everything in assembler, avoiding the deathly slow FPU, and doing trig by table lookup instead of blindly-repeated recalculation gave a speed improvement of over 100x... plenty fast enough for real-time on a 4.77 MHz PX/XT when coupled with fast assembly language display graphics.
Over the years, I have occasionally looked at touted "fast" FFT routines, but they appeared to all suffer from the same issues. (Like assuming that you were only going to perform one FFT, instead of pre-computing factors and swap tables so you don't have to repeat it for every FFT.)
Now, however, I suspect that much faster FPUs (and SSE, and GPUs) have made it actually better to recompute everything, rather than incur cache misses while trying to access precomputed tables. It seems that technology has finally vindicated the old "don't worry, someday the system will be fast enough" !
Best regards,
Bob Masta DAQARTA v6.00 Data AcQuisition And Real-Time Analysis
formatting link
Scope, Spectrum, Spectrogram, Sound Level Meter Frequency Counter, FREE Signal Generator Pitch Track, Pitch-to-MIDI Science with your sound card!
Thanks. I suppose in the design of radiotelescope arrays, folks keep an eye on the condition number of the transformation, at least at some level, and you know the Jacobian analytically, which is also a big help. Stuff like missing samples in a time series can give rise to all sorts of funnies if you try just interpolating.
Cheers
Phil Hobbs
--
Dr Philip C D Hobbs
Principal
ElectroOptical Innovations
55 Orchard Rd
Briarcliff Manor NY 10510
845-480-2058
email: hobbs (atsign) electrooptical (period) net
http://electrooptical.net
The history is quite interesting in that the interferometer arrays became more complex and completely beyond analytical treatment as computers became more powerful.
The first Earth rotation aperture synthesis telescope the One Mile Telescope at Cambridge had 2 fixed dishes and one movable on a very accurately surveyed E-W line. I think 64 baselines at uniform spacing (32 days observing) and 32 days moving the scope configuration.
The next generation 5km telescope was about 5 degrees of perfect E-W alignment (thanks to Dr Beeching closing the Cambridge-Bedford railway) and consisted of 4 fixed and 4 movable. 128 uniformly spaced baselines in 8 days observing and 8 days moving the scope around.
Big disadvantage of these E-W arrays is that their beamshape is lousy near the celestial equator and you have to observe for 12 hours continuous to get a decent image of the target. But on the plus side the mathematics for gridding the data is easier.
The VLA moved to a Y shape with 27 antennae I forget how many fixed and how many movable and a power law antenna spacing which allows more or less instant snapshots. When commissioned at the opening party it did some quick give away 5 minute maps of famous objects for VIPs to take away. It require massively more computing power to do this. The beamshape is pretty hairy and must be deconvolved by either CLEAN or MEM to provide an image or map suitable for astrophysics.
If you are interesting in the details a decent talk is online at
formatting link
It loses something by not having a narrator but if you are already familiar with the basics it should make reasonable sense. It shows what the U-V baseline coverage, beamshape, raw and final maps look like.
International VLBI you are pretty much stuck with where the biggest dishes are physically located on the planet and have to make do with what is one offer. The beamshape is horrible and it can only realistically be used on very bright compact targets. Cunning calibration techinques using triples of baselines extract good observables despite the unknown phase errors over each telescope. The latter method has been applied in the optical at near IR wavelengths.
Indeed. If you have missing data the only reliable way of modelling it is to pose the question what is the most probable model of the thing being measured that would yield the observed data to within the measurement noise. Attempts at direct inverses with missing data can end up dominated by the misleading zeroes or guessed missing data.
Astronomers do have some advantaged here. The sky brightness is known to be everywhere positive so negative regions have to be artefacts. The same cannot be said for electrical signals...
ElectronDepot website is not affiliated with any of the manufacturers or service providers discussed here.
All logos and trade names are the property of their respective owners.