PSD of a time series

Hi everyone,

I have a current control loop that drives a PWM and an Half-Bridge to
drive a BLDC motor phase. I measure current in the phase through a shunt
and an ADC fed by an OPA. The sampled data is fed back in the loop.

I'd like to calculate the PSD of the error between setpoint and measured
current at the input of the controller.

Assuming that I have my model working and that I have my time series for
the error signal, what would be the correct procedure to estimate the PSD?

Considering that the requirement is foreseen for a stable input (input
signal fixed), what is the best way to extract the PSD out of the model?
Should I maintain a fixed input? How do I decide the windowing? I know
that my requirement is in a specific frequency range and expressed in
mV^2/kHz.

A colleague of mine has performed the calculation by doing and FFT of
the signal and doing the pointwise product of the signal with the conjugate:

<m-code>
[spectrum,ffd] = fourierT(inSignals(:,1),sampling);
resBW=ffd(2)-ffd(1) %bandwidth of frequency bin
powerSpectrum=spectrum.*conj(spectrum)/resBW; %noise power spectrum
</m-code>

but the result is 3 order of magnitude *better* than the requirement,
which leads me to thinking that either we missed something or we've
overdesigned this thing by far!

Any idea/suggestion?

Al

-- 
A: Because it messes up the order in which people normally read text.
Q: Why is top-posting such a bad thing?
A: Top-posting.
Q: What is the most annoying thing on usenet and in e-mail?

On Fri, 03 Jun 2016 18:50:46 +0200, alb wrote:

> Hi everyone,
> 
> I have a current control loop that drives a PWM and an Half-Bridge to
> drive a BLDC motor phase. I measure current in the phase through a shunt
> and an ADC fed by an OPA. The sampled data is fed back in the loop.
> 
> I'd like to calculate the PSD of the error between setpoint and measured
> current at the input of the controller.
> 
> Assuming that I have my model working and that I have my time series for
> the error signal, what would be the correct procedure to estimate the
> PSD?
> 
> Considering that the requirement is foreseen for a stable input (input
> signal fixed), what is the best way to extract the PSD out of the model?
> Should I maintain a fixed input? How do I decide the windowing? I know
> that my requirement is in a specific frequency range and expressed in
> mV^2/kHz.
> 
> A colleague of mine has performed the calculation by doing and FFT of
> the signal and doing the pointwise product of the signal with the
> conjugate:
> 
> <m-code>
> [spectrum,ffd] = fourierT(inSignals(:,1),sampling);
> resBW=ffd(2)-ffd(1) %bandwidth of frequency bin
> powerSpectrum=spectrum.*conj(spectrum)/resBW; %noise power spectrum
> </m-code>
> 
> but the result is 3 order of magnitude *better* than the requirement,
> which leads me to thinking that either we missed something or we've
> overdesigned this thing by far!
> 
> Any idea/suggestion?
> 
> Al

I'm a bit of an amateur at PSD measurement, so take this all with a grain 
of salt.

Basically, when you take the Fourier Transform of a random sequence, you 
get a random sequence.  To get an estimate of the actual expected value 
of the PSD, you need to get the RMS average of a number of sample vectors 
(by taking a bunch of shorter FFTs and RMS-ing the results bin-by-bin).

That's all I know -- aside from knowing that I'd do energetic Googling 
and asking here.

-- 

Tim Wescott
Wescott Design Services
http://www.wescottdesign.com

I'm looking for work -- see my website!

You mention that you want the PSD of the error between the Pwm input and the current output for a DC input. You also mention that the result should be in mv^2 /kHz. So there is a long list of scale factors that influence this result. 

The value of the sense resistor
The vref of the adc
The gain of the pwm amp

Are you accounting for all these? And since your answer must be in terms of the psd of a voltage, what does that voltage correspond to? The voltage across the sense resistor?

Bob

The PSD is defined as the cosine transform of the autocorrelation
function, and this is often the simplest way to compute it.

Windowing enters into the picture twice: you can (and usually should)
window each segment of data before computing the autocorrelation;
and you can (and usually should) window the autocorrelation result
before taking its transform.

The length of the first window is mostly determined by how often you
need to compute the PSD (e.g. every 100 msec, every second, etc.).
It may also be a function of the coherence time of your signal, that
is, the interval over which the signal's second-order statistics 
are expected to be stationary.

The length of the second window, which is shorter than the first,
is determined by the required frequency resolution of the result,
and its length is on the order of the inverse of this resolution.

The type of the second window has significant impact on the nature
of the results.  It must have sufficient dynamic range for the 
application -- in this context, good "dynamic range", or stopband 
rejection, controls the extent to which energy at unrelated frequencies 
might spill into the frequencies of interest.  Popular windows include
the (standard) Hamming window, and higher-order Hamming windows that
exhibit more dynamic range.  But there are plenty of other windows,
Hann, Gaussian, etc.  Usually there would have to be something 
particular about your application to point you away from choosing
a Hamming window.

Hope this helps.

Steve

Hi Bob,

On 04.06.2016 13:45, radams2000@gmail.com wrote:
> You mention that you want the PSD of the error between the Pwm input
> and the current output for a DC input. You also mention that the
> result should be in mv^2 /kHz. So there is a long list of scale
> factors that influence this result.

I have a model of the loop, with all the related scale factors,
therefore my time series has already all these scale factor accounted
for. What is not taken into account is that the loop input and measured
current are expressed in Ampere, not Volts, or actually they are simply
a digital value that represent current. How should I change that to
Volts? Maybe it would be a better option to provide power spectrum in dB.

> Are you accounting for all these? And since your answer must be in
> terms of the psd of a voltage, what does that voltage correspond to?
> The voltage across the sense resistor?

That is a good question, since my setpoint is expressed in terms of
current, I suppose that my converted current will have conversion factor
that reflects the digitized quantity as current.

-- 
A: Because it messes up the order in which people normally read text.
Q: Why is top-posting such a bad thing?
A: Top-posting.
Q: What is the most annoying thing on usenet and in e-mail?

Hi Steve,

On 04.06.2016 14:47, Steve Pope wrote:
> The PSD is defined as the cosine transform of the autocorrelation
> function, and this is often the simplest way to compute it.
> 
> Windowing enters into the picture twice: you can (and usually should)
> window each segment of data before computing the autocorrelation;

Why should I use windowing for the autocorrelation?

> and you can (and usually should) window the autocorrelation result
> before taking its transform.

And why should I window the results for the transform?

> The length of the first window is mostly determined by how often you
> need to compute the PSD (e.g. every 100 msec, every second, etc.).

This is a measurement to assess the error distribution of our loop and
to verify that we meet the requirements, therefore is by nature
'offline' (no need to repeat it every n seconds). I assume then that the
windowing can be extended to the entire data sample. But how can I judge
how many points I need to include in my data sample?

> It may also be a function of the coherence time of your signal, that
> is, the interval over which the signal's second-order statistics 
> are expected to be stationary.

Uhm, I'm kinda lost. The PSD should be calculated in stationary
conditions, therefore the setpoint is set to a constant value 'C' and
the error w.r.t. the applied current is measured. I suspect I can let
the system 'stabilize' for a given amount of time (possibly several
times the 'settling time' of my transfer function) and then perform the
measurement.

> The length of the second window, which is shorter than the first,
> is determined by the required frequency resolution of the result,
> and its length is on the order of the inverse of this resolution.

That means that if I only have the range I wouldn't be able to assess
the PSD correctly, is that right? What if I blindly take the entire set
of sampled data again? Shouldn't I get all the information included in
my data sample?

> The type of the second window has significant impact on the nature
> of the results.  It must have sufficient dynamic range for the 
> application -- in this context, good "dynamic range", or stopband 
> rejection, controls the extent to which energy at unrelated frequencies 
> might spill into the frequencies of interest.  

Not sure I've grasped this. If my target requirement is a one sided PSD
of x mV^2/kHz over a frequency range between [0, 1.5 kHz], do you mean
that depending on the windowing I risk to have high frequencies
components 'spilling' into my [0, 1.5 kHz] bandwidth polluting my results?

> Popular windows include
> the (standard) Hamming window, and higher-order Hamming windows that
> exhibit more dynamic range.  But there are plenty of other windows,
> Hann, Gaussian, etc.  Usually there would have to be something 
> particular about your application to point you away from choosing
> a Hamming window.

I'm not really 100% confident in using a window without knowing if it is
the right method to use. What are those 'particular things' that may
drive me away from using Hamming windows? Any pointer?

-- 
A: Because it messes up the order in which people normally read text.
Q: Why is top-posting such a bad thing?
A: Top-posting.
Q: What is the most annoying thing on usenet and in e-mail?

On Sun, 05 Jun 2016 14:36:20 +0200, alb wrote:

> Hi Bob,
> 
> On 04.06.2016 13:45, radams2000@gmail.com wrote:
>> You mention that you want the PSD of the error between the Pwm input
>> and the current output for a DC input. You also mention that the result
>> should be in mv^2 /kHz. So there is a long list of scale factors that
>> influence this result.
> 
> I have a model of the loop, with all the related scale factors,
> therefore my time series has already all these scale factor accounted
> for. What is not taken into account is that the loop input and measured
> current are expressed in Ampere, not Volts, or actually they are simply
> a digital value that represent current. How should I change that to
> Volts? Maybe it would be a better option to provide power spectrum in
> dB.
> 
>> Are you accounting for all these? And since your answer must be in
>> terms of the psd of a voltage, what does that voltage correspond to?
>> The voltage across the sense resistor?
> 
> That is a good question, since my setpoint is expressed in terms of
> current, I suppose that my converted current will have conversion factor
> that reflects the digitized quantity as current.

If you have a quantity that's in amps, it would seem that you'd want a 
"power" spectrum in units of amps.

-- 

Tim Wescott
Wescott Design Services
http://www.wescottdesign.com

I'm looking for work -- see my website!

I'm curious about what you are looking for. If this is a PWM control loop where the input is DC and the system has already settled, then all you should see is a DC component (which should disappear if you subtract input from output) and a bunch of harmonics of the PWM carrier frequency. I assume your adc is anti-aliased so that you don't see the PWN components. So then ideally you should see nothing. 

Bob

On Sat, 4 Jun 2016 12:47:56 +0000 (UTC), spope33@speedymail.org (Steve
Pope) wrote:

>The PSD is defined as the cosine transform of the autocorrelation
>function, and this is often the simplest way to compute it.

That's one "definition", but there are many.   It is a pretty complex
topic.

A couple simple resources that delve into some of the nuances and
tradeoffs of different methods:

https://en.wikipedia.org/wiki/Spectral_density_estimation

And this has a bit of discussion around some of the tradeoffs of
methodologies:

http://dsp.stackexchange.com/questions/2096/why-so-many-methods-of-computing-psd

Many of the tradeoffs are application dependent, and I'm not sure I
understand the OPs app enough to add much more guidance.

Hi Bob,

On 05.06.2016 18:31, radams2000@gmail.com wrote:
> I'm curious about what you are looking for. If this is a PWM control
> loop where the input is DC and the system has already settled, then
> all you should see is a DC component (which should disappear if you
> subtract input from output) and a bunch of harmonics of the PWM
> carrier frequency. 

we're looking after noise. The loop has many sources of noise, from
current reading circuit noise, to ADC quantization, to H-bridge supply
noise. The aim of the exercise is to contain the noise in the system
within certain limits, that I'm not sure are correctly expressed in
mV^2/kHz at this point.

> I assume your adc is anti-aliased so that you
> don't see the PWN components.

The anti-aliasing cut-off is traded off with the phase shift introduced
in the loop, too much of it and I'd be risking instability. The high
frequency components of the PWM are filtered-out by an LC filter to
maintain common/differential mode emissions low.

> So then ideally you should see nothing.

The problem is that the world as we know it is never ideal ;-). My task
is to assess this distribution with some kind of analysis before
building the thing (no trials and errors allowed in a space program).
Now I have my models for the noise levels at each point in the loop,
estimated for the worst case and assumed to be independent. The customer
expressed his requirement as the one sided PSD and I'm assuming that is
one way to express noise figures, isn't it right?

Now that I have my model, I do not understand is what stimulus should I
inject? a constant one? and the for how many points should I measure the
output? All these points are kind of fuzzy today and I think I lack the
necessary background to find an answer to all this.

-- 
A: Because it messes up the order in which people normally read text.
Q: Why is top-posting such a bad thing?
A: Top-posting.
Q: What is the most annoying thing on usenet and in e-mail?