Groups | Search | Server Info | Keyboard shortcuts | Login | Register [http] [https] [nntp] [nntps]

Groups > comp.lang.python > #45557 > unrolled thread

Harmonic distortion of a input signal

Back to article view | Back to comp.lang.python

Contents

  Harmonic distortion of a input signal Anti Log <antilogeffects@gmail.com> - 2013-05-19 03:52 -0700
    Re: Harmonic distortion of a input signal Chris Angelico <rosuav@gmail.com> - 2013-05-20 01:34 +1000
      Re: Harmonic distortion of a input signal killybeard91@gmail.com - 2013-05-19 15:25 -0700
        Re: Harmonic distortion of a input signal Oscar Benjamin <oscar.j.benjamin@gmail.com> - 2013-05-19 23:49 +0100
        Re: Harmonic distortion of a input signal Terry Jan Reedy <tjreedy@udel.edu> - 2013-05-19 19:19 -0400
    Re: Harmonic distortion of a input signal killybeard91@gmail.com - 2013-05-19 15:59 -0700
    Re: Harmonic distortion of a input signal killybeard91@gmail.com - 2013-05-19 16:03 -0700
    Re: Harmonic distortion of a input signal killybeard91@gmail.com - 2013-05-19 16:36 -0700
      Re: Harmonic distortion of a input signal Gregory Ewing <greg.ewing@canterbury.ac.nz> - 2013-05-20 13:09 +1200
      Re: Harmonic distortion of a input signal Dave Angel <davea@davea.name> - 2013-05-19 21:11 -0400
        Re: Harmonic distortion of a input signal jmfauth <wxjmfauth@gmail.com> - 2013-05-20 10:23 -0700
          Re: Harmonic distortion of a input signal Christian Gollwitzer <auriocus@gmx.de> - 2013-05-20 19:50 +0200
            Re: Harmonic distortion of a input signal Christian Gollwitzer <auriocus@gmx.de> - 2013-05-20 19:56 +0200
              Re: Harmonic distortion of a input signal jmfauth <wxjmfauth@gmail.com> - 2013-05-23 04:44 -0700
          Re: Harmonic distortion of a input signal Oscar Benjamin <oscar.j.benjamin@gmail.com> - 2013-05-21 15:58 +0100

#45557 — Harmonic distortion of a input signal

Hi,
I have a task to calculate total distortion of a harmonics, of a signal that i imported from oscilloscope as numpy array. I had no problem drawing its spectrum, and time domain graph, but cant seem to find any functions that calculate TDH.
Any help? 
Best regards

[toc] | [next] | [standalone]

#45562

On Sun, May 19, 2013 at 8:52 PM, Anti Log <antilogeffects@gmail.com> wrote:
> total distortion of a harmonics

I selected this part of your post, right-clicked, and Chrome offered
to do a Google search for those words. And, surprise surprise, the
first hit is a page that appears to have the mathematics behind it.
Now, I don't know how much you trust Google and Wikipedia, but I'm
sure you can confirm the maths in some other way.

My guess is that there's no function in numpy to do what you're
asking... but it shouldn't be too hard to render the formula/e given
into Python code. Python's pretty expressive when it comes to algebra.
:)

ChrisA

[toc] | [prev] | [next] | [standalone]

#45568

How can i at least find a peek in FFT spectrum of a square wave ? 
From there i could easily build formula. Sorry for bothering but i am new to Python.

[toc] | [prev] | [next] | [standalone]

#45569

On 19 May 2013 23:25,  <killybeard91@gmail.com> wrote:
> How can i at least find a peek in FFT spectrum of a square wave ?
> From there i could easily build formula. Sorry for bothering but i am new to Python.

Are you the same person who posted the original question?

You probably want to use numpy for this. I'm not sure if I understand
your question but here goes:

First import numpy (you may need to install this first):

>>> import numpy as np

Create a square wave signal:

>>> x = np.zeros(50)
>>> x[:25] = -1
>>> x[25:] = +1
>>> x
array([-1., -1., -1., -1., -1., -1., -1., -1., -1., -1., -1., -1., -1.,
       -1., -1., -1., -1., -1., -1., -1., -1., -1., -1., -1., -1.,  1.,
        1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,
        1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.])

Compute the magnitude spectrum:

>>> spect = abs(np.fft.fft(x)[:25])
>>> spect
array([  0.        ,  31.85194222,   0.        ,  10.67342282,
         0.        ,   6.47213595,   0.        ,   4.69726931,
         0.        ,   3.73254943,   0.        ,   3.13762901,
         0.        ,   2.7436023 ,   0.        ,   2.47213595,
         0.        ,   2.28230601,   0.        ,   2.15105461,
         0.        ,   2.06487174,   0.        ,   2.01589594,   0.        ])

Find the index of the maximum element:

>>> np.argmax(spect)
1

So the peak is the lowest non-zero frequency component of the DFT. In
Hz this corresponds to a frequency of 1/T where T is the duration of
the signal.


Oscar

[toc] | [prev] | [next] | [standalone]

#45572

On 5/19/2013 6:49 PM, Oscar Benjamin wrote:

>>>> import numpy as np
>
> Create a square wave signal:
>
>>>> x = np.zeros(50)
>>>> x[:25] = -1
>>>> x[25:] = +1
>>>> x
> array([-1., -1., -1., -1., -1., -1., -1., -1., -1., -1., -1., -1., -1.,
>         -1., -1., -1., -1., -1., -1., -1., -1., -1., -1., -1., -1.,  1.,
>          1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,
>          1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.])
>
> Compute the magnitude spectrum:
>
>>>> spect = abs(np.fft.fft(x)[:25])
>>>> spect
> array([  0.        ,  31.85194222,   0.        ,  10.67342282,
>           0.        ,   6.47213595,   0.        ,   4.69726931,
>           0.        ,   3.73254943,   0.        ,   3.13762901,
>           0.        ,   2.7436023 ,   0.        ,   2.47213595,
>           0.        ,   2.28230601,   0.        ,   2.15105461,
>           0.        ,   2.06487174,   0.        ,   2.01589594,   0.        ])
>
> Find the index of the maximum element:
>
>>>> np.argmax(spect)
> 1
>
> So the peak is the lowest non-zero frequency component of the DFT. In
> Hz this corresponds to a frequency of 1/T where T is the duration of
> the signal.

While you were answering a specific question, I think the above is a 
nice tutorial example, because it is more meaningful than arbitrary 
operations applied to random data.

[toc] | [prev] | [next] | [standalone]

#45570

Yes, sorry logged from another account.
Would that work on a numpy array ? Because this signal was imported from oscilloscope as a numpy array.
Best regards,

[toc] | [prev] | [next] | [standalone]

#45571

Got it working, thanks alot :)

[toc] | [prev] | [next] | [standalone]

#45573

One more question. Function np.argmax returns max of non-complex numbers ? 
Because FFT array of my signal is complex.

[toc] | [prev] | [next] | [standalone]

#45575

killybeard91@gmail.com wrote:
> One more question. Function np.argmax returns max of non-complex numbers ? 
> Because FFT array of my signal is complex.

You'll want the magnitudes of the complex numbers.
Actually the squares of the magnitudes (assuming the
data from the oscilloscope represents voltages),
because you're after a ratio of powers.

-- 
Greg

[toc] | [prev] | [next] | [standalone]

#45576

On 05/19/2013 07:36 PM, killybeard91@gmail.com wrote:
> One more question. Function np.argmax returns max of non-complex numbers ?
> Because FFT array of my signal is complex.
>

It'd be easier to track the thread if you actually replied to the 
message you're responding to, and also if you included some context. 
But I'll paste the latter in here:

Terry Reedy said:
 > Compute the magnitude spectrum:

 >>> spect = abs(np.fft.fft(x)[:25])
 >>> spect
 > array([  0.        ,  31.85194222,   0.        ,  10.67342282,
 >          0.        ,   6.47213595,   0.        ,   4.69726931,
 >          0.        ,   3.73254943,   0.        ,   3.13762901,
 >          0.        ,   2.7436023 ,   0.        ,   2.47213595,
 >          0.        ,   2.28230601,   0.        ,   2.15105461,
 >          0.        ,   2.06487174,   0.        ,   2.01589594,
 > 0.        ])

 > Find the index of the maximum element:

 >>> np.argmax(spect)
 > 1


Notice that argmax's argument is the result of an abs() call.  It's got 
real numbers representing the magnitude of the various complex numbers.

-- 
DaveA

[toc] | [prev] | [next] | [standalone]

#45625

Non sense.

The discrete fft algorithm is valid only if the number of data
points you transform does correspond to a power of 2 (2**n).

Keywords to the problem: apodization, zero filling, convolution
product, ...

eg. http://en.wikipedia.org/wiki/Convolution

jmf

[toc] | [prev] | [next] | [standalone]

#45626

Am 20.05.13 19:23, schrieb jmfauth:
> Non sense.

Dito.

> The discrete fft algorithm is valid only if the number of data
> points you transform does correspond to a power of 2 (2**n).

Where did you get this? The DFT is defined for any integer point number 
the same way.

Just if you want to get it fast, you need to worry about the length. For 
powers of two, there is the classic Cooley-Tukey. But there do exist FFT 
algorithms for any other length. For example, there is the Winograd 
transform for a set of small numbers, there is "mixed-radix" to reduce 
any length which can be factored, and there is finally Bluestein which 
works for any size, even for a prime. All of the aforementioned 
algorithms are O(log n) and are implemented in typical FFT packages. All 
of them should result (up to rounding differences) in the same thing as 
the naive DFT sum. Therefore, today

> Keywords to the problem: apodization, zero filling, convolution
> product, ...

Not for a periodic signal of integer length.

> eg. http://en.wikipedia.org/wiki/Convolution

How long do you read this group?

	Christian

[toc] | [prev] | [next] | [standalone]

#45627

Oops, I thought we were posting to comp.dsp. Nevertheless, I think 
numpy.fft does mixed-radix (can't check it now)

Am 20.05.13 19:50, schrieb Christian Gollwitzer:
> Am 20.05.13 19:23, schrieb jmfauth:
>> Non sense.
>
> Dito.
>
>> The discrete fft algorithm is valid only if the number of data
>> points you transform does correspond to a power of 2 (2**n).
>
> Where did you get this? The DFT is defined for any integer point number
> the same way.
>
> Just if you want to get it fast, you need to worry about the length. For
> powers of two, there is the classic Cooley-Tukey. But there do exist FFT
> algorithms for any other length. For example, there is the Winograd
> transform for a set of small numbers, there is "mixed-radix" to reduce
> any length which can be factored, and there is finally Bluestein which
> works for any size, even for a prime. All of the aforementioned
> algorithms are O(log n) and are implemented in typical FFT packages. All
> of them should result (up to rounding differences) in the same thing as
> the naive DFT sum. Therefore, today
>
>> Keywords to the problem: apodization, zero filling, convolution
>> product, ...
>
> Not for a periodic signal of integer length.
>
>> eg. http://en.wikipedia.org/wiki/Convolution
>
> How long do you read this group?
>
>      Christian
>

[toc] | [prev] | [next] | [standalone]

#45802

On 20 mai, 19:56, Christian Gollwitzer <aurio...@gmx.de> wrote:
> Oops, I thought we were posting to comp.dsp. Nevertheless, I think
> numpy.fft does mixed-radix (can't check it now)
>
> Am 20.05.13 19:50, schrieb Christian Gollwitzer:
>
>
>
>
>
>
>
> > Am 20.05.13 19:23, schrieb jmfauth:
> >> Non sense.
>
> > Dito.
>
> >> The discrete fft algorithm is valid only if the number of data
> >> points you transform does correspond to a power of 2 (2**n).
>
> > Where did you get this? The DFT is defined for any integer point number
> > the same way.
>
> > Just if you want to get it fast, you need to worry about the length. For
> > powers of two, there is the classic Cooley-Tukey. But there do exist FFT
> > algorithms for any other length. For example, there is the Winograd
> > transform for a set of small numbers, there is "mixed-radix" to reduce
> > any length which can be factored, and there is finally Bluestein which
> > works for any size, even for a prime. All of the aforementioned
> > algorithms are O(log n) and are implemented in typical FFT packages. All
> > of them should result (up to rounding differences) in the same thing as
> > the naive DFT sum. Therefore, today
>
> >> Keywords to the problem: apodization, zero filling, convolution
> >> product, ...
>
> > Not for a periodic signal of integer length.
>
> >> eg.http://en.wikipedia.org/wiki/Convolution
>
> > How long do you read this group?
>
> >      Christian

------

Forget what I wrote.
I'm understanding what I wanted to say, it is badly
formulated.

jmf

[toc] | [prev] | [next] | [standalone]

#45669

On 20 May 2013 18:23, jmfauth <wxjmfauth@gmail.com> wrote:
> Non sense.
>
> The discrete fft algorithm is valid only if the number of data
> points you transform does correspond to a power of 2 (2**n).

As with many of your comments about Python's unicode implementation
you are confusing performance with validity. The DFT is defined and is
a valid invertible map (barring roundoff) for complex vectors of any
integer length. It is also a valid method for understanding the
frequency content of periodic signals. The fastest FFT algorithms are
for vectors whose length is a power of 2 but the other algorithms
produce equally *valid* DFT results.

In the example I posted the computation of the DFT using numpy.fft.fft
was (as far as I could tell) instantaneous. I could use timeit to
discover exactly how many microseconds it took but why when I already
have the results I wanted?

> Keywords to the problem: apodization, zero filling, convolution
> product, ...
>
> eg. http://en.wikipedia.org/wiki/Convolution

These points are not relevant to the example given.


Oscar

[toc] | [prev] | [standalone]

Back to top | Article view | comp.lang.python

csiph-web