On 10/16/2011 5:33 PM, mnentwig wrote:
>>> A good general rule:
> ..
>
> Engineering time is very expensive those days. CPU time comes for free in
> small amounts.
> A 10M elements FFT occupies one of the four cores in my desktop PC for ~500
> ms.
> If you see any valid reasons to use more engineering time to improve on
> that number, fine.
> I don't.
>
> If it were about selling 100 million gadgets that implement the algorithm,
> we could talk about it. But here, it's about getting the job done in the
> lab, I suppose.
>
> Besides, if I get too clever in my implementation, it becomes inflexible.
> Why should I create my own maintenance nightmare, when I can keep it
> simple.

No argument. I was just telling you how I think. I spent too many years 
doing tasks with an 1802 that "experts" claimed required at least an 80186.

Jerry
-- 
Engineering is the art of making what you want from things you can get.

>> A good general rule: 
..

Engineering time is very expensive those days. CPU time comes for free in
small amounts.
A 10M elements FFT occupies one of the four cores in my desktop PC for ~500
ms. 
If you see any valid reasons to use more engineering time to improve on
that number, fine. 
I don't.

If it were about selling 100 million gadgets that implement the algorithm,
we could talk about it. But here, it's about getting the job done in the
lab, I suppose.

Besides, if I get too clever in my implementation, it becomes inflexible.
Why should I create my own maintenance nightmare, when I can keep it
simple.

On 10/15/2011 3:21 AM, mnentwig wrote:
> Hello,

   ...

> What I don't understand are some of the responses to the OP's question in
> general.
> What on earth is wrong with using FFT - I need 10 lines of code, do we
> expect a more elegant solution still?

I tend to think in terms of machine cycles rather than lines of code. 
After all, with the right stuff behind the curtain, the code can be 
boiled down to
   run(f,a)
or some such.

> The accuracy of the frequency response is perfect, that's hard to beat. The
> only catch is, I get a line spectrum, from the periodicity of the signal.
> But the measurement duration is finite in any case =>  increase the size as
> needed.

Prtfection is overkill.

> What's the problem with generating Gaussian noise in the frequency domain
> using a Gaussian amplitude and uniform [0..2pi[ phase per frequency bin?
> This will work just fine (use conj() on negative frequencies, as was said
> earlier, to get a real-valued signal in the time domain).

The problem? too many chances to go wrong.

A good general rule: If the code is for something mechanical, you can 
almost always do the job with an 8-bit processor running at less than 10 
MHz. Before I get swamped, remember that I wrote 'can', not 'should'.

Jerry
-- 
Engineering is the art of making what you want from things you can get.

>> The AMPLITUDE should be Rayleigh distributed
thanks, I agree.

On Oct 15, 3:21�am, "mnentwig"
<markus.nentwig@n_o_s_p_a_m.renesasmobile.com> wrote:
<<<material snipped>>>
>
> What's the problem with generating Gaussian noise in the frequency domain
> using a Gaussian amplitude and uniform [0..2pi[ phase per frequency bin?
> This will work just fine (use conj() on negative frequencies, as was said
> earlier, to get a real-valued signal in the time domain).


If you use a Gaussian AMPLITUDE and uniform phase in the
frequency domain and taking the iDFT, then you are relying
on the central limit theorem to give you Gaussian noise in
the time domain.  (The AMPLITUDE should be Rayleigh
distributed).  What you should be doing is putting in complex
Gaussian noise in each bin (zero-mean equal-variance
independent Gaussian RVs in real and imaginary parts,
conjugates in the negative frequencies).  This is not too hard
in general.  The Box-Mulller method that is commonly used
returns TWO independent standard Gaussian random variables
on each call, exactly what is needed except for  a scaling factor.

Dilip Sarwate

>> in the frequency domain using a Gaussian amplitude and uniform [0..2pi[
phase per frequency bin?7

better: ... uniform [0..pi[ phase, since the Gaussian amplitude can be
negative. 

-Markus

PS: To clarify, referring to Tim's first post:

The key is:
"Any time domain transients of the processing simply wrap around into the
future / past cycle." 

As long as I don't care about time domain information, this is not a
problem at all (as said, tested and tried in the lab, accurate to
instrument precision). 
A noise signal does NOT contain any timing information or the like, so
we're on the safe side.

Now If I'd try the same method with a _modulated_ signal, I need to be
aware that the infinitely long impulse response simply wraps around.
I observe one cycle of a periodic signal, with the decaying transients of
all past cycles and the rising transients of all future cycles overlapping.
The result is correct (time-domain convolution of a periodic signal with an
infinite impulse response). Whether it's what I want is another question,
but I can easily modify the frequency response (raised-cosine, for example,
comes to mind).

Hello,

>So what's the big deal?  It is known that a
>Gaussian random process passing through
>a linear filter emerges as a Gaussian process,
>and so the histogram of values of sig *should*
>look Gaussian, and by golly it does. ...

Thanks, that's exactly what I meant.

> so
>you could equally well have started by filling
>sig with 475,000 Gaussian random variables
>at one end, their complex conjugates at the
>other end, leaving 50,000 zeroes in the middle
>and taken the iDFT or ifft.  

I agree. 
In one way, I'd explain it intuitively as "filtering" (=> remove energy
beyond a cutoff frequency), in the other way as sample rate conversion. The
outcome is the same.

>> And I won't bring up
>the problems that Tim Wescott and Glenn
>Herrmansfeldt have already noted with this
>brick-wall filter business.

There is not necessarily a problem, if used in the OP's context (driving a
vibration exciter). 
I've been using this method for years in the lab with IQ modulation signals
for RF signal generators. Those signals have to be periodic anyway, since
the memory size on the instruments used to be limited. Therefore, FFT is a
perfect tool IMHO: Any time domain transients of the processing simply wrap
around into the future / past cycle. 

If I set a brick-wall spectrum on the generated signal, I see _exactly_ a
brick wall spectrum on spectrum analyzer or on the downconverted and
sampled signal back in Matlab. Just what I'd expect.

What I don't understand are some of the responses to the OP's question in
general. 
What on earth is wrong with using FFT - I need 10 lines of code, do we
expect a more elegant solution still? 
The accuracy of the frequency response is perfect, that's hard to beat. The
only catch is, I get a line spectrum, from the periodicity of the signal.
But the measurement duration is finite in any case => increase the size as
needed.

What's the problem with generating Gaussian noise in the frequency domain
using a Gaussian amplitude and uniform [0..2pi[ phase per frequency bin?
This will work just fine (use conj() on negative frequencies, as was said
earlier, to get a real-valued signal in the time domain).

-Markus

On Oct 13, 6:45�pm, "mnentwig"
<markus.nentwig@n_o_s_p_a_m.renesasmobile.com>
set down a program and remarked about its output:

> Looks gaussian enough to me.

I know just enough programming to make some sense
of what he wrote, so I am going to put down what I think
the program does and then make some additional
comments based on my understanding

> sig = randn(1, n);
Fills the array sig with the values returned by n = 10^6
calls to a standard Gaussian random variable generator
randn.  The histogram of 10^6 values will look like a
standard Gaussian pdf

> sig = fft(sig);
Compute in-place DFT of the vector sig of length 10^6

> sig (1+n2:end-n2+1) = 0;
Zero out the middle 50,000 bin, effectively passing
sig through a brick-wall filter

> sig = ifft(sig);
Compute in-place iDFT of filtered signal

> sig = real(sig); % chop off roundoff error
Discard imaginary components since filtered
output should be a real-valued signal

> hist(sig, 100); hold on;
> x = -3:0.01:3;
> plot(x, 36000*exp(-(2.2*x) .^ 2), 'r');
Show the histogram of sig and overlay a
Gaussian pdf curve over it to show that
the histogram pretty much follows the same
shape.

====================================
So what's the big deal?  It is known that a
Gaussian random process passing through
a linear filter emerges as a Gaussian process,
and so the histogram of values of sig *should*
look Gaussian, and by golly it does.  So what?
As Jerry Avins has already pointed out, a lot
of work has been done for little effect; sig
began as a Gaussian vector anyway.  Oh, you
mean that sig is now a low-pass Gaussian
process as opposed to the broadband Gaussian
process that you began with?  Well, why didn't
you say so?  But please note also that
>sig =  fft(sig) is also a Gaussian vector, so
you could equally well have started by filling
sig with 475,000 Gaussian random variables
at one end, their complex conjugates at the
other end, leaving 50,000 zeroes in the middle
and taken the iDFT or ifft.  And I won't bring up
the problems that Tim Wescott and Glenn
Herrmansfeldt have already noted with this
brick-wall filter business.

Dilip Sarwate

On 10/14/2011 2:09 PM, mnentwig wrote:
> yes, but isn't it supposed to be Gaussian?
> What the FFT does is simple oversampling.

It isn't exactly Gaussian, but very close. (No value exceeds .5) See the 
central limit theorem.

> If the discrete-time samples have a Gaussian PDF, shouldn't the same apply
> for the continuous-time waveform which is reconstructed by an ideal lowpass
> filter (or approximated by the FFT with arbitrary accuracy by using more
> zero padding)?

Do you think an FFT takes less computation than adding 12 random numbers 
spit out by a linear congruential generator? How much accuracy do you 
expect a shaker table to need?

Jerry
--                                           the least expensive
Engineering is the art of making what you want from things you can get.
                                                    ^