frequency components of a signal buried in a noisy time domain signal
Hi all, A common use of Fourier transforms is to find the frequency components of a signal buried in a noisy time domain signal. I found a Matlab template at http://www.mathworks.com/access/helpdesk/help/techdoc/ref/fft.shtml Matlab has a function nextpow2(L) Is there a similar build-in function in numpy/scipy ? I tried to convert the m-file into a pythonic form. What is needed to obtain a similar figure of the single-sided amplitude spectrum using numpy/scipy/matplotlib ? from numpy import sin, linspace, pi from numpy.random import randn from pylab import plot, show, title, xlabel, ylabel from scipy.fft import fft Fs = 1000. # Sampling frequency T = 1./Fs # Sample time L = 1000 # length of signal t = arange(0,L)*T x = 0.7*sin(2*pi*50*t)+sin(2*pi*120*t) y = x + 2*randn(len(t)) plot(Fs*t[:50],y[:50]) title('Signal corrupted with zero-mean random noise') xlabel('Time (milliseconds)') # # # #NFFT = 2^nextpow2(L); # Next power of 2 from length of y Y = fft(y,NFFT)/L f = Fs/2*linspace(0,1,NFFT/2+1) plot(f,2*abs(Y[:NFFT/2+1])) title('Single-sided amplitude spectrum of y(t)') xlabel('Frequency (Hz)') ylabel('|Y(f)|') show() What can be done in case of nonequispaced data ? http://dx.doi.org/10.1137/0914081 Thanks in advance Nils
Nils, Yes, FFT can be used to visualize the frequency content of a signal buried in the noise, especially if it is is narrowband, even though more advanced spectral analysis is required for advanced applications. I cannot help you with non-uniformly spaced samples (some sort of interpolation comes to mind), but next power of 2 should be trivial. Without checking for the argument type, you can implement the function like this: def nextpow2(n): m_f = np.log2(n) m_i = np.ceil(m_f) return 2**m_i Hope it helps, Ivo On 26 February 2010 14:05, Nils Wagner <nwagner@iam.uni-stuttgart.de> wrote:
Hi all,
A common use of Fourier transforms is to find the frequency components of a signal buried in a noisy time domain signal.
I found a Matlab template at
http://www.mathworks.com/access/helpdesk/help/techdoc/ref/fft.shtml
Matlab has a function
nextpow2(L)
Is there a similar build-in function in numpy/scipy ?
I tried to convert the m-file into a pythonic form. What is needed to obtain a similar figure of the single-sided amplitude spectrum using numpy/scipy/matplotlib ?
from numpy import sin, linspace, pi from numpy.random import randn from pylab import plot, show, title, xlabel, ylabel from scipy.fft import fft
Fs = 1000. # Sampling frequency T = 1./Fs # Sample time L = 1000 # length of signal t = arange(0,L)*T x = 0.7*sin(2*pi*50*t)+sin(2*pi*120*t) y = x + 2*randn(len(t)) plot(Fs*t[:50],y[:50]) title('Signal corrupted with zero-mean random noise') xlabel('Time (milliseconds)') # # # #NFFT = 2^nextpow2(L); # Next power of 2 from length of y
Y = fft(y,NFFT)/L f = Fs/2*linspace(0,1,NFFT/2+1) plot(f,2*abs(Y[:NFFT/2+1])) title('Single-sided amplitude spectrum of y(t)') xlabel('Frequency (Hz)') ylabel('|Y(f)|') show()
What can be done in case of nonequispaced data ?
http://dx.doi.org/10.1137/0914081
Thanks in advance
Nils _______________________________________________ SciPy-User mailing list SciPy-User@scipy.org http://mail.scipy.org/mailman/listinfo/scipy-user
On Fri, 26 Feb 2010 15:56:29 -0500 Ivo Maljevic <ivo.maljevic@gmail.com> wrote:
Nils, Yes, FFT can be used to visualize the frequency content of a signal buried in the noise, especially if it is is narrowband, even though more advanced spectral analysis is required for advanced applications.
Thank you very much for your reply. Are you aware of free software for advanced applications ? What could be done in case of broadband noise ? Nils
I cannot help you with non-uniformly spaced samples (some sort of interpolation comes to mind), but next power of 2 should be trivial. Without checking for the argument type, you can implement the function like this:
def nextpow2(n): m_f = np.log2(n) m_i = np.ceil(m_f) return 2**m_i
Hope it helps, Ivo
On 26 February 2010 14:05, Nils Wagner <nwagner@iam.uni-stuttgart.de> wrote:
Hi all,
A common use of Fourier transforms is to find the frequency components of a signal buried in a noisy time domain signal.
I found a Matlab template at
http://www.mathworks.com/access/helpdesk/help/techdoc/ref/fft.shtml
Matlab has a function
nextpow2(L)
Is there a similar build-in function in numpy/scipy ?
I tried to convert the m-file into a pythonic form. What is needed to obtain a similar figure of the single-sided amplitude spectrum using numpy/scipy/matplotlib ?
from numpy import sin, linspace, pi from numpy.random import randn from pylab import plot, show, title, xlabel, ylabel from scipy.fft import fft
Fs = 1000. # Sampling frequency T = 1./Fs # Sample time L = 1000 # length of signal t = arange(0,L)*T x = 0.7*sin(2*pi*50*t)+sin(2*pi*120*t) y = x + 2*randn(len(t)) plot(Fs*t[:50],y[:50]) title('Signal corrupted with zero-mean random noise') xlabel('Time (milliseconds)') # # # #NFFT = 2^nextpow2(L); # Next power of 2 from length of y
Y = fft(y,NFFT)/L f = Fs/2*linspace(0,1,NFFT/2+1) plot(f,2*abs(Y[:NFFT/2+1])) title('Single-sided amplitude spectrum of y(t)') xlabel('Frequency (Hz)') ylabel('|Y(f)|') show()
What can be done in case of nonequispaced data ?
http://dx.doi.org/10.1137/0914081
Thanks in advance
Nils _______________________________________________ SciPy-User mailing list SciPy-User@scipy.org http://mail.scipy.org/mailman/listinfo/scipy-user
Hi, Looking at a periodic signal buried in noise is a well-studied problem, with many techniques for attacking it. You really need to be a little more specific about what you want to do. For example, is your input signal really a sinusoid, or does it have harmonic content? Are you trying to detect a weak periodic signal or are you trying to extract the features of a strong periodic signal? Is your signal exactly periodic, does it have some (deterministic or random) wander, or are you looking for the power spectrum of a broadband signal? If your input data are non-uniformly sampled, everything becomes more difficult (and computationally expensive), but there are solutions (e.g. the Lomb-Scargle periodogram). Anne On 26 February 2010 16:24, Nils Wagner <nwagner@iam.uni-stuttgart.de> wrote:
On Fri, 26 Feb 2010 15:56:29 -0500 Ivo Maljevic <ivo.maljevic@gmail.com> wrote:
Nils, Yes, FFT can be used to visualize the frequency content of a signal buried in the noise, especially if it is is narrowband, even though more advanced spectral analysis is required for advanced applications.
Thank you very much for your reply.
Are you aware of free software for advanced applications ? What could be done in case of broadband noise ?
Nils
I cannot help you with non-uniformly spaced samples (some sort of interpolation comes to mind), but next power of 2 should be trivial. Without checking for the argument type, you can implement the function like this:
def nextpow2(n): m_f = np.log2(n) m_i = np.ceil(m_f) return 2**m_i
Hope it helps, Ivo
On 26 February 2010 14:05, Nils Wagner <nwagner@iam.uni-stuttgart.de> wrote:
Hi all,
A common use of Fourier transforms is to find the frequency components of a signal buried in a noisy time domain signal.
I found a Matlab template at
http://www.mathworks.com/access/helpdesk/help/techdoc/ref/fft.shtml
Matlab has a function
nextpow2(L)
Is there a similar build-in function in numpy/scipy ?
I tried to convert the m-file into a pythonic form. What is needed to obtain a similar figure of the single-sided amplitude spectrum using numpy/scipy/matplotlib ?
from numpy import sin, linspace, pi from numpy.random import randn from pylab import plot, show, title, xlabel, ylabel from scipy.fft import fft
Fs = 1000. # Sampling frequency T = 1./Fs # Sample time L = 1000 # length of signal t = arange(0,L)*T x = 0.7*sin(2*pi*50*t)+sin(2*pi*120*t) y = x + 2*randn(len(t)) plot(Fs*t[:50],y[:50]) title('Signal corrupted with zero-mean random noise') xlabel('Time (milliseconds)') # # # #NFFT = 2^nextpow2(L); # Next power of 2 from length of y
Y = fft(y,NFFT)/L f = Fs/2*linspace(0,1,NFFT/2+1) plot(f,2*abs(Y[:NFFT/2+1])) title('Single-sided amplitude spectrum of y(t)') xlabel('Frequency (Hz)') ylabel('|Y(f)|') show()
What can be done in case of nonequispaced data ?
http://dx.doi.org/10.1137/0914081
Thanks in advance
Nils _______________________________________________ SciPy-User mailing list SciPy-User@scipy.org http://mail.scipy.org/mailman/listinfo/scipy-user
_______________________________________________ SciPy-User mailing list SciPy-User@scipy.org http://mail.scipy.org/mailman/listinfo/scipy-user
On Fri, 26 Feb 2010 16:32:01 -0500 Anne Archibald <peridot.faceted@gmail.com> wrote:
Hi,
Looking at a periodic signal buried in noise is a well-studied problem, with many techniques for attacking it. You really need to be a little more specific about what you want to do. For example, is your input signal really a sinusoid, or does it have harmonic content? Are you trying to detect a weak periodic signal or are you trying to extract the features of a strong periodic signal? Is your signal exactly periodic, does it have some (deterministic or random) wander, or are you looking for the power spectrum of a broadband signal?
If your input data are non-uniformly sampled, everything becomes more difficult (and computationally expensive), but there are solutions (e.g. the Lomb-Scargle periodogram).
Anne
Hi Anne, Thank you very much for your hints ! BTW, a BSD licensed code for the Lomb-Scargle periodogram is available at http://www.mathworks.com/matlabcentral/fileexchange/993-lombscargle-m http://www.mathworks.com/matlabcentral/fileexchange/20004-lomb-lomb-scargle-... I am newbie to signal processing. Is there a good introduction that you can recommend ? There are so many books on signal processing. It should cover engineering applications. What makes a signal weak/strong periodic ? The signals come from real-life application (pressure / acceleration data). Do I need a filter before I apply FFT ? What would you do if you know nothing about the origin of the signal ? How can I distinguish between deterministic and random wander ? Cheers, Nils
Nils, I will not answer your questions to Anne then (I could, but you asked her :) ), but I'll give you another advice, given that you are new to signal processing. Whatever you do, please do not use the scipy's snr function call to calculate the SNR in signal processing. I'll illustrate the problem with the 'statistical SNR' in signal processing with 3 closely related and simple examples. 1. Infinitely long sinusoidal signal with additive white noise Let x(t)=s(t)+n(t)=A sin(2\pi f t) + n(t) Power of the signal component can be found as P_s=1/(2T) \int_T^T s^2(t) dt = A^2/2, and let the power of noise be P_n=var(n) = 1. If, we chose A=100, then intuitively, we know that the SNR is large, and indeed: SNR = P_s/ P_N = A^2/2 = 5000 On the other hand, if you use statistical definition (used in optical astronomy, and scipy and implemented in SCIPY), you would conclude quite the opposite: SNR = \mu_s / signa_n = mean(s) / signa_n= 0/1 = 0, and that would be the case regardless of what the value for A is. 2. Discrete time signal with finite duration Let's take this same signal and sample it in a such way that we have 100 cycles of a sinusoid:
f=200 t=np.linspace(0,0.5,2000) s=100*np.sin(2*np.pi*f*t)
Using sampled signal definitions for powers, you get: P_s = np.mean(s**2) = 4997.5 and you can calculate the power of the noise either as:
np.var(n) = 0.97726608241925972
or the same way as the power of the signal:
np.mean(n**2) = 0.977977713
The results will slightly differ, and that has to do with biased/unbiased estimator. The SNR=4997.5/0.9772660=5113.756 which is pretty close to analog domain case with infinitely long signal. Using the scipy approach, you would however conclude the following: SNR=mean(s)/sigma_n=-5.6495516860933771e-15/0.9885676923808808= -5.7148860210947362e-15 (which comes close to zero). 3. Same example, but let's take fractional number of cycles.
t=np.linspace(0,0.5+np.pi/800,2000+int(5*np.pi)) # this is 100 cycles + 1/4 of a cycle s=100*np.sin(2*np.pi*f*t)
P_s = np.mean(s**2) SNR=P_s/np.var(n) SNR 5117.8701015651177 # very close to the two results above
But, using the other definition:
SNR=np.mean(s)/np.std(n) SNR 0.098933308385744489
As you can see, your SNR would be a function of the frequency of the signal (that is, it would matter how you truncate the cycles) and it wouldn't matter much what is the amplitude of the signal. Just my 2 cents, Ivo On 27 February 2010 04:41, Nils Wagner <nwagner@iam.uni-stuttgart.de> wrote:
On Fri, 26 Feb 2010 16:32:01 -0500 Anne Archibald <peridot.faceted@gmail.com> wrote:
Hi,
Looking at a periodic signal buried in noise is a well-studied problem, with many techniques for attacking it. You really need to be a little more specific about what you want to do. For example, is your input signal really a sinusoid, or does it have harmonic content? Are you trying to detect a weak periodic signal or are you trying to extract the features of a strong periodic signal? Is your signal exactly periodic, does it have some (deterministic or random) wander, or are you looking for the power spectrum of a broadband signal?
If your input data are non-uniformly sampled, everything becomes more difficult (and computationally expensive), but there are solutions (e.g. the Lomb-Scargle periodogram).
Anne
Hi Anne,
Thank you very much for your hints !
BTW, a BSD licensed code for the Lomb-Scargle periodogram is available at http://www.mathworks.com/matlabcentral/fileexchange/993-lombscargle-m
http://www.mathworks.com/matlabcentral/fileexchange/20004-lomb-lomb-scargle-...
I am newbie to signal processing. Is there a good introduction that you can recommend ? There are so many books on signal processing. It should cover engineering applications.
What makes a signal weak/strong periodic ?
The signals come from real-life application (pressure / acceleration data). Do I need a filter before I apply FFT ?
What would you do if you know nothing about the origin of the signal ?
How can I distinguish between deterministic and random wander ?
Cheers, Nils _______________________________________________ SciPy-User mailing list SciPy-User@scipy.org http://mail.scipy.org/mailman/listinfo/scipy-user
On Sat, Feb 27, 2010 at 6:41 PM, Nils Wagner <nwagner@iam.uni-stuttgart.de> wrote:
On Fri, 26 Feb 2010 16:32:01 -0500 Anne Archibald <peridot.faceted@gmail.com> wrote:
Hi,
Looking at a periodic signal buried in noise is a well-studied problem, with many techniques for attacking it. You really need to be a little more specific about what you want to do. For example, is your input signal really a sinusoid, or does it have harmonic content? Are you trying to detect a weak periodic signal or are you trying to extract the features of a strong periodic signal? Is your signal exactly periodic, does it have some (deterministic or random) wander, or are you looking for the power spectrum of a broadband signal?
If your input data are non-uniformly sampled, everything becomes more difficult (and computationally expensive), but there are solutions (e.g. the Lomb-Scargle periodogram).
Anne
Hi Anne,
Thank you very much for your hints !
BTW, a BSD licensed code for the Lomb-Scargle periodogram is available at http://www.mathworks.com/matlabcentral/fileexchange/993-lombscargle-m http://www.mathworks.com/matlabcentral/fileexchange/20004-lomb-lomb-scargle-...
I am newbie to signal processing. Is there a good introduction that you can recommend ?
It depends on what you are looking for, and the depth and time you are willing to spend on it. A "down to earth" book is freely available here: http://www.dspguide.com If you are interested in manipulating signals with noise, statistical signal processing is where to look at, but I am not sure there is any good book which does not start with theory there.
There are so many books on signal processing. It should cover engineering applications.
What makes a signal weak/strong periodic ?
The signals come from real-life application (pressure / acceleration data). Do I need a filter before I apply FFT ?
What would you do if you know nothing about the origin of the signal ?
Well, if you know really nothing, there is nothing you can do :) For spectral estimation, there are broadly two kind of methods: parametric and non parametric. Unfortunately, scipy itself does not have much here - I started the talkbox scikits to implement at least what is available in matlab and more for spectral estimation, but I am not sure I will have time to work on it much anymore. Basic methods, like AR modelling, are quite simple to implement by yourself, though. David
On Sat, Feb 27, 2010 at 7:44 AM, David Cournapeau <cournape@gmail.com> wrote:
On Sat, Feb 27, 2010 at 6:41 PM, Nils Wagner <nwagner@iam.uni-stuttgart.de> wrote:
On Fri, 26 Feb 2010 16:32:01 -0500 Anne Archibald <peridot.faceted@gmail.com> wrote:
Hi,
Looking at a periodic signal buried in noise is a well-studied problem, with many techniques for attacking it. You really need to be a little more specific about what you want to do. For example, is your input signal really a sinusoid, or does it have harmonic content? Are you trying to detect a weak periodic signal or are you trying to extract the features of a strong periodic signal? Is your signal exactly periodic, does it have some (deterministic or random) wander, or are you looking for the power spectrum of a broadband signal?
If your input data are non-uniformly sampled, everything becomes more difficult (and computationally expensive), but there are solutions (e.g. the Lomb-Scargle periodogram).
Anne
Hi Anne,
Thank you very much for your hints !
BTW, a BSD licensed code for the Lomb-Scargle periodogram is available at http://www.mathworks.com/matlabcentral/fileexchange/993-lombscargle-m http://www.mathworks.com/matlabcentral/fileexchange/20004-lomb-lomb-scargle-...
I am newbie to signal processing. Is there a good introduction that you can recommend ?
It depends on what you are looking for, and the depth and time you are willing to spend on it. A "down to earth" book is freely available here:
If you are interested in manipulating signals with noise, statistical signal processing is where to look at, but I am not sure there is any good book which does not start with theory there.
There are so many books on signal processing. It should cover engineering applications.
What makes a signal weak/strong periodic ?
The signals come from real-life application (pressure / acceleration data). Do I need a filter before I apply FFT ?
What would you do if you know nothing about the origin of the signal ?
Well, if you know really nothing, there is nothing you can do :) For spectral estimation, there are broadly two kind of methods: parametric and non parametric. Unfortunately, scipy itself does not have much here - I started the talkbox scikits to implement at least what is available in matlab and more for spectral estimation, but I am not sure I will have time to work on it much anymore.
Basic methods, like AR modelling, are quite simple to implement by yourself, though.
nitime also has many functions for frequency domain analysis and there are also a few functions in matplotlib Josef
David _______________________________________________ SciPy-User mailing list SciPy-User@scipy.org http://mail.scipy.org/mailman/listinfo/scipy-user
On 27 February 2010 04:41, Nils Wagner <nwagner@iam.uni-stuttgart.de> wrote:
On Fri, 26 Feb 2010 16:32:01 -0500
I am newbie to signal processing. Is there a good introduction that you can recommend ? There are so many books on signal processing. It should cover engineering applications.
Hmm. I don't have a good example in mind. For an overview you could try Numerical Recipes (but keep in mind that there are free python implementations of almost all the algorithms they describe, and there are often better algorithms out there). Mostly I have used signal processing in a scientific context, often somewhat different.
What makes a signal weak/strong periodic ?
By an exactly periodic signal I mean something like sin(f*t). Such a signal produces a very narrow peak of a characteristic shape in an FFT, and so can be recognized even in the presence of quite strong noise. If your signal is only quasi-periodic - perhaps the frequency is a slowly-varying function of time - you'll have a much broader peak, which will be much lower for the same input signal amplitude, and hence more difficult to distinguish from noise. If your signal is only quasi-periodic, or comes and goes in the data, you may want to do a series of FFTs on short, overlapping pieces of the data, so you can look at time evolution of the signal's spectral properties. By weak versus strong, I just meant: is your signal so weak that there's some question you'll be able to detect it at all, or is it strong enough that you can try to extract more subtle properties?
The signals come from real-life application (pressure / acceleration data). Do I need a filter before I apply FFT ?
The short answer is "probably not", or at least, not in the usual sense of "filter". You will probably want to extend your data set so that it is a power of two, just to make the FFT faster. Extending your data set in time has the effect of interpolating the FFT, so it may well be worth extending the data set further, to double or four times its original length. This will broaden your peaks and affect the statistics of the noise. (And, incidentally, if you're going to do this, you should extend the array either by repeating the input values or by filling the extra space with the mean of the input array to avoid a spurious bump at low frequencies.) Once you've got the FFT, you should know what different kinds of signal look like. A pure sinusoid that's at one of the FFT frequencies - which happens when there's exactly an integer number of turns in your data - will prodduce a spike at exactly one of your FFT frequencies. A pure sinusoid that's between two FFT frequencies will produce a sin(x)/x response in nearby bins; in particular, a sinusoid midway between two FFT frequencies will produce two high-power bins that are substantially less high than the single bin would be if the same signal were exactly at an FFT frequency; this is one reason extending the data to interpolate the FFT can be a good idea. There are also tricks that can be done once you've found a peak, looking at the precise shape of the peak to determine whether the periodic signal fills the whole of your data and if not where it is. If you have a more broad-band signal, perhaps because it is only approximately periodic, then you'll have an extended region of extra power in the FFT. This can be a little tricky to recognize, if it's not obvious from inspection; one thing to try is smoothing the power spectrum (the squared amplitude of the FFT) so that local increases show more prominently.
What would you do if you know nothing about the origin of the signal ?
This rather depends on what I was looking for. But if I had a time series in which I suspected some sort of periodicity, I'd start by plotting the power spectrum and looking for peaks that stood out well above the local average. I'd also try plotting a smoothed power spectrum. Also worth trying would be taking a series of overlapping FFTs each covering part of the data; I might see a signal in some of these individual FFTs, which would tell me where it was in the data, or I might average all the power spectra (which is roughly equivalent to smoothing a long-term FFT) looking for a broad-band signal. If I expected a periodic signal that had harmonic content - a series of pulses rather than a pure sinusoid, for example - I'd want to look at the total power in several harmonics at once, adding the power at f, 2f, 3f, 4f, .... If I found an exactly periodic signal, I could investigate further by "folding" the data, that is, adding all the data points based on the phase at which they were taken to build up a "profile" of the signal as a function of pulse phase.
How can I distinguish between deterministic and random wander ?
Without some a priori knowledge, it's not easy, though overlapping short-time FFTs might help. Here I'd produce a grayscale plot of power versus frequency and time (i.e. which piece the short FFT was applied to). If you have a nice bright signal, you might see it at the same frequency all the time, or you might see it drift linearly up or down, or it might vary sinusoidally in frequency over time, or it might wander randomly, or it might just always be broad. If you know exactly how it should be varying, say frequency varying linearly with time (e.g. Doppler shift from a uniform acceleration) or sinusoidally, there are special techniques to remove or detect that modulation. Some of these techniques are done in the time domain; others by looking at the shape of the peak (if any) in the FFT. I should say that this is a big enough and complex enough field that my answers reflect only my experience, which is skewed towards working with pulsar data, in which the data tends to include very weak but exactly periodic signals. People used to working with, say, speech data will have very different ideas on what sort of tools are useful. Anne
What makes a signal weak/strong periodic ?
By an exactly periodic signal I mean something like sin(f*t). Such a signal produces a very narrow peak of a characteristic shape in an FFT, and so can be recognized even in the presence of quite strong noise. If your signal is only quasi-periodic - perhaps the frequency is a slowly-varying function of time - you'll have a much broader peak, which will be much lower for the same input signal amplitude, and hence more difficult to distinguish from noise. If your signal is only quasi-periodic, or comes and goes in the data, you may want to do a series of FFTs on short, overlapping pieces of the data, so you can look at time evolution of the signal's spectral properties.
In my opinion this is not quite so. Periodic signal, as you rightly pointed out, is a sinusoidal signal. Quasi-periodic signal behaves like a periodic one, even though it does not satisfy the periodic condition x(t) = x (t+To), where To is the period. Best known examples are when you add two sinusoidal signals with frequencies that are not a fractional integer of each other. For example: sin(2pi f t)+sin(2pi^2 f t). You would still see a "spike" in the frequency domain, but quasi-periodicity definitely does not relate to low frequency. Ivo
On 27 February 2010 08:37, Ivo Maljevic <ivo.maljevic@gmail.com> wrote:
What makes a signal weak/strong periodic ?
By an exactly periodic signal I mean something like sin(f*t). Such a signal produces a very narrow peak of a characteristic shape in an FFT, and so can be recognized even in the presence of quite strong noise. If your signal is only quasi-periodic - perhaps the frequency is a slowly-varying function of time - you'll have a much broader peak, which will be much lower for the same input signal amplitude, and hence more difficult to distinguish from noise. If your signal is only quasi-periodic, or comes and goes in the data, you may want to do a series of FFTs on short, overlapping pieces of the data, so you can look at time evolution of the signal's spectral properties.
In my opinion this is not quite so. Periodic signal, as you rightly pointed out, is a sinusoidal signal. Quasi-periodic signal behaves like a periodic one, even though it does not satisfy the periodic condition x(t) = x (t+To), where To is the period. Best known examples are when you add two sinusoidal signals with frequencies that are not a fractional integer of each other. For example: sin(2pi f t)+sin(2pi^2 f t). You would still see a "spike" in the frequency domain, but quasi-periodicity definitely does not relate to low frequency.
I guess we have a minor conflict of definitions. Of course the signal you describe is not periodic, but, at least in terms of Fourier transforms, you can treat the terms of a sum as completely independent signals. So I would view this as simply two exactly periodic signals, which you can indeed pull out as two separate peaks in the Fourier transform. What I had in mind in referring to quasi-periodic signals was something you might get out of an oscillator with poor frequency stability, or a not-very-good resonant cavity, or a sinusoid with some modulation: a signal localized in frequency but not exactly periodic. This will have a broader peak in the Fourier transform, but will still be interesting in many contexts. Anne
I guess we have a minor conflict of definitions. Of course the signal you describe is not periodic, but, at least in terms of Fourier transforms, you can treat the terms of a sum as completely independent signals. So I would view this as simply two exactly periodic signals, which you can indeed pull out as two separate peaks in the Fourier transform.
You are right. I just checked, and what I was referring to is not known as * quasi-periodic* signal, it is known as *almost-periodic*, and it plays a role when you want to discuss power vs. energy signals. Therefore, I take back my comment, as I am not familiar with the quasi-periodic signal definition. I am certainly familiar with phase noise / frequency error behaviour of oscillators, but that doesn't matter. Since you didn't recommend any book, I'll recommend a couple of books to Nils: 1. http://www.amazon.com/Discrete-Time-Signal-Processing-2nd-Prentice-Hall/dp/0... 2. http://www.amazon.com/Digital-Signal-Processing-Bookware-Companion/dp/049507... Ivo
Maybe you can try sigview (http://www.sigview.com/), if that is what you need, or you can invest some time in learning more about spectral analysis, and then write your own software that will be tied to your particular application / requirement. I cannot give you a short answer to your question regarding the noise bandwidth, but I can give you a few pointers. 1. When you have sampled signal which consists of desired signal + noise, usually it helps to know what is the bandwidth of the desired signal. In most applications that is the case, unless you are analyzing a completely unknown signal (military applications, SETI, etc). Then you can filter the signal to remove as much noise as possible. 2. In real world applications, there will always be some wideband noise (i.e., thermal noise), and you might also have some narrowband noise/interference. White noise is the most common type of wideband noise that you would encounter, even though it may be "coloured" noise if you are analyzing a signal at the output of nonideal Rx filter (for example in radio systems). 3. As long as power spectral density of the signal is higher then the power spectral density of the noise (at least in some portions of the spectrum), simple FFT or periodogram (do a google search) approach will suffice. Cheers, Ivo On 26 February 2010 16:24, Nils Wagner <nwagner@iam.uni-stuttgart.de> wrote:
On Fri, 26 Feb 2010 15:56:29 -0500 Ivo Maljevic <ivo.maljevic@gmail.com> wrote:
Nils, Yes, FFT can be used to visualize the frequency content of a signal buried in the noise, especially if it is is narrowband, even though more advanced spectral analysis is required for advanced applications.
Thank you very much for your reply.
Are you aware of free software for advanced applications ? What could be done in case of broadband noise ?
Nils
I cannot help you with non-uniformly spaced samples (some sort of interpolation comes to mind), but next power of 2 should be trivial. Without checking for the argument type, you can implement the function like this:
def nextpow2(n): m_f = np.log2(n) m_i = np.ceil(m_f) return 2**m_i
Hope it helps, Ivo
On 26 February 2010 14:05, Nils Wagner <nwagner@iam.uni-stuttgart.de> wrote:
Hi all,
A common use of Fourier transforms is to find the frequency components of a signal buried in a noisy time domain signal.
I found a Matlab template at
http://www.mathworks.com/access/helpdesk/help/techdoc/ref/fft.shtml
Matlab has a function
nextpow2(L)
Is there a similar build-in function in numpy/scipy ?
I tried to convert the m-file into a pythonic form. What is needed to obtain a similar figure of the single-sided amplitude spectrum using numpy/scipy/matplotlib ?
from numpy import sin, linspace, pi from numpy.random import randn from pylab import plot, show, title, xlabel, ylabel from scipy.fft import fft
Fs = 1000. # Sampling frequency T = 1./Fs # Sample time L = 1000 # length of signal t = arange(0,L)*T x = 0.7*sin(2*pi*50*t)+sin(2*pi*120*t) y = x + 2*randn(len(t)) plot(Fs*t[:50],y[:50]) title('Signal corrupted with zero-mean random noise') xlabel('Time (milliseconds)') # # # #NFFT = 2^nextpow2(L); # Next power of 2 from length of y
Y = fft(y,NFFT)/L f = Fs/2*linspace(0,1,NFFT/2+1) plot(f,2*abs(Y[:NFFT/2+1])) title('Single-sided amplitude spectrum of y(t)') xlabel('Frequency (Hz)') ylabel('|Y(f)|') show()
What can be done in case of nonequispaced data ?
http://dx.doi.org/10.1137/0914081
Thanks in advance
Nils _______________________________________________ SciPy-User mailing list SciPy-User@scipy.org http://mail.scipy.org/mailman/listinfo/scipy-user
_______________________________________________ SciPy-User mailing list SciPy-User@scipy.org http://mail.scipy.org/mailman/listinfo/scipy-user
Nils, I modified a little bit your code. Try this version, maybe you will find it helpful. I wouldn't bother much with the next power of 2 unless you are going to implement something in real time - FFT calculation with non-power of 2 on a PC is just fine. Also, I would consider periodograms only if the signal changes over time, and you want to be able to track these changes over time. Otherwise, taking the whole signal and doing one FFT may be OK. Take a look at the code below (blue) , which is a modified version of your code (and read the comments). Ivo import numpy as np import matplotlib.pyplot as plt from scipy.fftpack import fft, fftshift def signal_spectrum(x, Fs): n = len(x) f = np.linspace(-Fs/2, Fs/2, n) X=abs(fftshift(fft(x))) # move the upper spectrum to negative freqs return [f, X] pi = np.pi Fs = 1000. # Sampling frequency T = 1./Fs # Sample time L = 1000 # length of signal t = np.arange(0,L)*T x = 0.7*np.sin(2*pi*50*t)+np.sin(2*pi*120*t) y = x + 2*np.random.randn(len(t)) plt.figure() plt.plot(Fs*t[:50],y[:50]) plt.title('Signal corrupted with zero-mean random noise') plt.xlabel('Time (milliseconds)') NFFT = 2**np.ceil(np.log2(L)) # only if you really want this plt.figure() Y = fft(y,NFFT)/L f = Fs/2*np.linspace(0,1,NFFT/2+1) plt.plot(f,2*abs(Y[:NFFT/2+1])) plt.title('Single-sided amplitude spectrum of y(t)') plt.xlabel('Frequency (Hz)') plt.ylabel('|Y(f)|') # Another way of calculating and plotting the spectrum f,Yd = signal_spectrum(y, Fs) plt.figure() # you want to add this to be able to see multiple plots plt.plot(f,Yd) plt.title('Spectrum of the the signal') plt.xlabel('Frequency (Hz)') plt.ylabel('|Y(f)|') plt.grid(True) plt.show() On 26 February 2010 14:05, Nils Wagner <nwagner@iam.uni-stuttgart.de> wrote:
Hi all,
A common use of Fourier transforms is to find the frequency components of a signal buried in a noisy time domain signal.
I found a Matlab template at
http://www.mathworks.com/access/helpdesk/help/techdoc/ref/fft.shtml
Matlab has a function
nextpow2(L)
Is there a similar build-in function in numpy/scipy ?
I tried to convert the m-file into a pythonic form. What is needed to obtain a similar figure of the single-sided amplitude spectrum using numpy/scipy/matplotlib ?
from numpy import sin, linspace, pi from numpy.random import randn from pylab import plot, show, title, xlabel, ylabel from scipy.fft import fft
Fs = 1000. # Sampling frequency T = 1./Fs # Sample time L = 1000 # length of signal t = arange(0,L)*T x = 0.7*sin(2*pi*50*t)+sin(2*pi*120*t) y = x + 2*randn(len(t)) plot(Fs*t[:50],y[:50]) title('Signal corrupted with zero-mean random noise') xlabel('Time (milliseconds)') # # # #NFFT = 2^nextpow2(L); # Next power of 2 from length of y
Y = fft(y,NFFT)/L f = Fs/2*linspace(0,1,NFFT/2+1) plot(f,2*abs(Y[:NFFT/2+1])) title('Single-sided amplitude spectrum of y(t)') xlabel('Frequency (Hz)') ylabel('|Y(f)|') show()
What can be done in case of nonequispaced data ?
http://dx.doi.org/10.1137/0914081
Thanks in advance
Nils _______________________________________________ SciPy-User mailing list SciPy-User@scipy.org http://mail.scipy.org/mailman/listinfo/scipy-user
On Fri, 26 Feb 2010 22:03:08 -0500 Ivo Maljevic <ivo.maljevic@gmail.com> wrote:
Nils, I modified a little bit your code. Try this version, maybe you will find it helpful. I wouldn't bother much with the next power of 2 unless you are going to implement something in real time - FFT calculation with non-power of 2 on a PC is just fine. Also, I would consider periodograms only if the signal changes over time, and you want to be able to track these changes over time. Otherwise, taking the whole signal and doing one FFT may be OK.
Take a look at the code below (blue) , which is a modified version of your code (and read the comments).
Ivo
import numpy as np import matplotlib.pyplot as plt from scipy.fftpack import fft, fftshift
def signal_spectrum(x, Fs): n = len(x) f = np.linspace(-Fs/2, Fs/2, n) X=abs(fftshift(fft(x))) # move the upper spectrum to negative freqs return [f, X]
pi = np.pi Fs = 1000. # Sampling frequency T = 1./Fs # Sample time L = 1000 # length of signal t = np.arange(0,L)*T x = 0.7*np.sin(2*pi*50*t)+np.sin(2*pi*120*t) y = x + 2*np.random.randn(len(t)) plt.figure() plt.plot(Fs*t[:50],y[:50]) plt.title('Signal corrupted with zero-mean random noise') plt.xlabel('Time (milliseconds)')
NFFT = 2**np.ceil(np.log2(L)) # only if you really want this
plt.figure() Y = fft(y,NFFT)/L f = Fs/2*np.linspace(0,1,NFFT/2+1) plt.plot(f,2*abs(Y[:NFFT/2+1])) plt.title('Single-sided amplitude spectrum of y(t)') plt.xlabel('Frequency (Hz)') plt.ylabel('|Y(f)|')
# Another way of calculating and plotting the spectrum f,Yd = signal_spectrum(y, Fs) plt.figure() # you want to add this to be able to see multiple plots plt.plot(f,Yd) plt.title('Spectrum of the the signal') plt.xlabel('Frequency (Hz)') plt.ylabel('|Y(f)|') plt.grid(True)
plt.show()
Hi Ivo, Thank you very much for your hints. My experience is, that spectrograms work well with synthetic signals. What is needed to use spectrograms for real-life signals ? I mean do I need a filter etc. ? Nils from scipy import linspace, vectorize, cos, pi from pylab import plot, show, imshow, subplot, specgram, xlabel, ylabel, savefig # # Example taken from http://en.wikipedia.org/wiki/Short-time_Fourier_transform # def x(t): " Synthetic signal " if t < 5: return cos(2*pi*10*t) if t >= 5. and t < 10: return cos(2*pi*25*t) if t >=10. and t< 15: return cos(2*pi*50*t) if t >=15. and t<= 20: return cos(2*pi*100*t) t = linspace(0.,20.,8001) # sampled at 400 Hz dt = t[1]-t[0] NFFT = 512 # the length of the windowing segments Fs = int(1.0/dt) # the sampling frequency x_vec = vectorize(x) signal = x_vec(t) ax1=subplot(211) plot(t,signal) subplot(212,sharex=ax1) Pxx, freqs, bins, im = specgram(signal, NFFT=NFFT, Fs=Fs, noverlap=5) xlabel(r'Time $t$ [s]') ylabel('Frequency $f$ [Hz]') savefig('spectrogram',dpi=60) show()
Hi Ivo,
Thank you very much for your hints. My experience is, that spectrograms work well with synthetic signals.
What is needed to use spectrograms for real-life signals ?
The same way the frequency content of your synthetic signal changes over time, it happens with real signal (e.g., voice, communication system's signal, etc). I mean do I need a filter etc. ?
There are several levels of filtering: 1. You want to filter a signal even before you sample it (e.g., in comm systems, you want to remove as much noise+interference) and focus only on the band that interests you, plus you want to satisfy Nyquist. 2. After you have sampled your signal (if it was analog to begin with, otherwise, skip step 1), you may do additional digital filtering, to further reduce the impact of noise/interference, but this is all application specific. 3. Depending on what you are doing, you might want to consider windowing techniques in your spectral analysis, if you want to suppress the sidelobes (fft uses rectangular window, which has high sidelobes - sinc function), or to better identify the peaks. It all comes down to trade-offs. In any case, to think about these problems, you need to go over signal processing material first. Ivo
Nils
On Fri, Feb 26, 2010 at 12:05 PM, Nils Wagner <nwagner@iam.uni-stuttgart.de>wrote:
Hi all,
A common use of Fourier transforms is to find the frequency components of a signal buried in a noisy time domain signal.
I found a Matlab template at
http://www.mathworks.com/access/helpdesk/help/techdoc/ref/fft.shtml
Matlab has a function
nextpow2(L)
I use searchsorted for that. In [5]: tab = 2**arange(30) In [6]: tab[tab.searchsorted(10)] Out[6]: 16 In [7]: tab[tab.searchsorted(81)] Out[7]: 128 Is there a similar build-in function in numpy/scipy ?
It would be handy sometimes. Where would it go? Chuck
On Sat, Feb 27, 2010 at 12:16 AM, Charles R Harris < charlesr.harris@gmail.com> wrote:
On Fri, Feb 26, 2010 at 12:05 PM, Nils Wagner < nwagner@iam.uni-stuttgart.de> wrote:
Hi all,
A common use of Fourier transforms is to find the frequency components of a signal buried in a noisy time domain signal.
I found a Matlab template at
http://www.mathworks.com/access/helpdesk/help/techdoc/ref/fft.shtml
Matlab has a function
nextpow2(L)
I use searchsorted for that.
In [5]: tab = 2**arange(30)
In [6]: tab[tab.searchsorted(10)] Out[6]: 16
In [7]: tab[tab.searchsorted(81)] Out[7]: 128
Is there a similar build-in function in numpy/scipy ?
It would be handy sometimes. Where would it go?
Why wouldn't it go in fftpack? DG
Chuck
_______________________________________________ SciPy-User mailing list SciPy-User@scipy.org http://mail.scipy.org/mailman/listinfo/scipy-user
David, Nice way of avoiding log2, but how do you determine the length of your tab array? Do you want to chose an arbitrary large array, or you still want to use log2 call based on an input argument to the function nexpow2()? In other words, how would you make a generic nextpow2() function? Ivo On 27 February 2010 03:22, David Goldsmith <d.l.goldsmith@gmail.com> wrote:
On Sat, Feb 27, 2010 at 12:16 AM, Charles R Harris < charlesr.harris@gmail.com> wrote:
I use searchsorted for that.
In [5]: tab = 2**arange(30)
In [6]: tab[tab.searchsorted(10)] Out[6]: 16
In [7]: tab[tab.searchsorted(81)] Out[7]: 128
Is there a similar build-in function in numpy/scipy ?
It would be handy sometimes. Where would it go?
Why wouldn't it go in fftpack?
DG
Chuck
_______________________________________________ SciPy-User mailing list SciPy-User@scipy.org http://mail.scipy.org/mailman/listinfo/scipy-user
_______________________________________________ SciPy-User mailing list SciPy-User@scipy.org http://mail.scipy.org/mailman/listinfo/scipy-user
On Sat, Feb 27, 2010 at 04:02, Ivo Maljevic <ivo.maljevic@gmail.com> wrote:
David, Nice way of avoiding log2, but how do you determine the length of your tab array?
The input arrays can only be so large. Even on 64-bit machines, the table need not have more than 64 entries. tab = 2 ** np.arange(np.iinfo(np.intp).bits) -- Robert Kern "I have come to believe that the whole world is an enigma, a harmless enigma that is made terrible by our own mad attempt to interpret it as though it had an underlying truth." -- Umberto Eco
On Sat, Feb 27, 2010 at 11:02 AM, Robert Kern <robert.kern@gmail.com> wrote:
On Sat, Feb 27, 2010 at 04:02, Ivo Maljevic <ivo.maljevic@gmail.com> wrote:
David, Nice way of avoiding log2, but how do you determine the length of your tab array?
The input arrays can only be so large. Even on 64-bit machines, the table need not have more than 64 entries.
tab = 2 ** np.arange(np.iinfo(np.intp).bits)
Or tab = 1 << arange(iinfo(intp).bits - 1) Chuck
On Sat, Feb 27, 2010 at 1:31 PM, Charles R Harris <charlesr.harris@gmail.com> wrote:
On Sat, Feb 27, 2010 at 11:02 AM, Robert Kern <robert.kern@gmail.com> wrote:
On Sat, Feb 27, 2010 at 04:02, Ivo Maljevic <ivo.maljevic@gmail.com> wrote:
David, Nice way of avoiding log2, but how do you determine the length of your tab array?
The input arrays can only be so large. Even on 64-bit machines, the table need not have more than 64 entries.
tab = 2 ** np.arange(np.iinfo(np.intp).bits)
Or
tab = 1 << arange(iinfo(intp).bits - 1)
I think I prefer a readable solution to this for calls that are not in an inner loop. scipy.signal.fftconvolve uses the same as Ivo's solution << means it's much larger than 1 ? Josef
Chuck
_______________________________________________ SciPy-User mailing list SciPy-User@scipy.org http://mail.scipy.org/mailman/listinfo/scipy-user
On Sat, Feb 27, 2010 at 12:41, <josef.pktd@gmail.com> wrote:
On Sat, Feb 27, 2010 at 1:31 PM, Charles R Harris <charlesr.harris@gmail.com> wrote:
On Sat, Feb 27, 2010 at 11:02 AM, Robert Kern <robert.kern@gmail.com> wrote:
On Sat, Feb 27, 2010 at 04:02, Ivo Maljevic <ivo.maljevic@gmail.com> wrote:
David, Nice way of avoiding log2, but how do you determine the length of your tab array?
The input arrays can only be so large. Even on 64-bit machines, the table need not have more than 64 entries.
tab = 2 ** np.arange(np.iinfo(np.intp).bits)
Or
tab = 1 << arange(iinfo(intp).bits - 1)
I think I prefer a readable solution to this for calls that are not in an inner loop.
scipy.signal.fftconvolve uses the same as Ivo's solution
<< means it's much larger than 1 ?
It's a bitshift operator. Basically "x << y" means "x * (2 ** y)" for integer arguments. -- Robert Kern "I have come to believe that the whole world is an enigma, a harmless enigma that is made terrible by our own mad attempt to interpret it as though it had an underlying truth." -- Umberto Eco
I am sorry guys, but the solution I've proposed, and the one that was proposed by the scipy camp would cause confusion for those who switch from Matlab, and I just checked that. If you are going to use the same name as the Matlab one, than this function should return the m such that 2**m >= abs(n), and should not return 2^m. Whatever approach you guys want to use (and I believe that especially given what the function should return), do not return 2**m, but just m instead. For example:
def nextpow2(n): ... return np.ceil(np.log2(n)) ... m=nextpow2(1000) m 10.0 2**m 1024.0
On 27 February 2010 13:45, Robert Kern <robert.kern@gmail.com> wrote:
On Sat, Feb 27, 2010 at 12:41, <josef.pktd@gmail.com> wrote:
On Sat, Feb 27, 2010 at 1:31 PM, Charles R Harris <charlesr.harris@gmail.com> wrote:
On Sat, Feb 27, 2010 at 11:02 AM, Robert Kern <robert.kern@gmail.com>
wrote:
On Sat, Feb 27, 2010 at 04:02, Ivo Maljevic <ivo.maljevic@gmail.com> wrote:
David, Nice way of avoiding log2, but how do you determine the length of
your
tab array?
The input arrays can only be so large. Even on 64-bit machines, the table need not have more than 64 entries.
tab = 2 ** np.arange(np.iinfo(np.intp).bits)
Or
tab = 1 << arange(iinfo(intp).bits - 1)
I think I prefer a readable solution to this for calls that are not in an inner loop.
scipy.signal.fftconvolve uses the same as Ivo's solution
<< means it's much larger than 1 ?
It's a bitshift operator. Basically "x << y" means "x * (2 ** y)" for integer arguments.
-- Robert Kern
"I have come to believe that the whole world is an enigma, a harmless enigma that is made terrible by our own mad attempt to interpret it as though it had an underlying truth." -- Umberto Eco _______________________________________________ SciPy-User mailing list SciPy-User@scipy.org http://mail.scipy.org/mailman/listinfo/scipy-user
On Sat, Feb 27, 2010 at 11:58 AM, Ivo Maljevic <ivo.maljevic@gmail.com>wrote:
I am sorry guys, but the solution I've proposed, and the one that was proposed by the scipy camp would cause confusion for those who switch from Matlab, and I just checked that. If you are going to use the same name as the Matlab one, than this function should return the m such that 2**m >= abs(n), and should not return 2^m.
Whatever approach you guys want to use (and I believe that especially given what the function should return), do not return 2**m, but just m instead. For example:
For the table approach, it is the index returned by searchsorted. For floats it is np.frexp(x)[1]. Chuck
Another suggestion, and this one also makes the parallel to Matlab. Make sure it works for vectors:
a=np.array([7,11,250]) nextpow2(a) array([ 3., 4., 8.])
On 27 February 2010 13:58, Ivo Maljevic <ivo.maljevic@gmail.com> wrote:
I am sorry guys, but the solution I've proposed, and the one that was proposed by the scipy camp would cause confusion for those who switch from Matlab, and I just checked that. If you are going to use the same name as the Matlab one, than this function should return the m such that 2**m >= abs(n), and should not return 2^m.
Whatever approach you guys want to use (and I believe that especially given what the function should return), do not return 2**m, but just m instead. For example:
def nextpow2(n): ... return np.ceil(np.log2(n)) ... m=nextpow2(1000) m 10.0 2**m 1024.0
On 27 February 2010 13:45, Robert Kern <robert.kern@gmail.com> wrote:
On Sat, Feb 27, 2010 at 12:41, <josef.pktd@gmail.com> wrote:
On Sat, Feb 27, 2010 at 1:31 PM, Charles R Harris <charlesr.harris@gmail.com> wrote:
On Sat, Feb 27, 2010 at 11:02 AM, Robert Kern <robert.kern@gmail.com>
wrote:
On Sat, Feb 27, 2010 at 04:02, Ivo Maljevic <ivo.maljevic@gmail.com> wrote:
David, Nice way of avoiding log2, but how do you determine the length of
your
tab array?
The input arrays can only be so large. Even on 64-bit machines, the table need not have more than 64 entries.
tab = 2 ** np.arange(np.iinfo(np.intp).bits)
Or
tab = 1 << arange(iinfo(intp).bits - 1)
I think I prefer a readable solution to this for calls that are not in an inner loop.
scipy.signal.fftconvolve uses the same as Ivo's solution
<< means it's much larger than 1 ?
It's a bitshift operator. Basically "x << y" means "x * (2 ** y)" for integer arguments.
-- Robert Kern
"I have come to believe that the whole world is an enigma, a harmless enigma that is made terrible by our own mad attempt to interpret it as though it had an underlying truth." -- Umberto Eco _______________________________________________ SciPy-User mailing list SciPy-User@scipy.org http://mail.scipy.org/mailman/listinfo/scipy-user
On Sat, Feb 27, 2010 at 12:03 PM, Ivo Maljevic <ivo.maljevic@gmail.com>wrote:
Another suggestion, and this one also makes the parallel to Matlab. Make sure it works for vectors:
a=np.array([7,11,250]) nextpow2(a) array([ 3., 4., 8.])
In [27]: tab.searchsorted([7,11,250]) Out[27]: array([3, 4, 8]) In [28]: np.frexp([7,11,250])[1] Out[28]: array([3, 4, 8], dtype=int32) <snip> Chuck
On Sat, Feb 27, 2010 at 2:07 PM, Charles R Harris <charlesr.harris@gmail.com> wrote:
On Sat, Feb 27, 2010 at 12:03 PM, Ivo Maljevic <ivo.maljevic@gmail.com> wrote:
Another suggestion, and this one also makes the parallel to Matlab. Make sure it works for vectors:
a=np.array([7,11,250]) nextpow2(a) array([ 3., 4., 8.])
In [27]: tab.searchsorted([7,11,250]) Out[27]: array([3, 4, 8])
In [28]: np.frexp([7,11,250])[1] Out[28]: array([3, 4, 8], dtype=int32)
It looks like frexp is doubling the array if n is already an integer power
np.frexp(2**np.arange(5)) (array([ 0.5, 0.5, 0.5, 0.5, 0.5]), array([1, 2, 3, 4, 5])) np.arange(5) array([0, 1, 2, 3, 4]) np.ceil(np.log2(2**np.arange(5))) array([ 0., 1., 2., 3., 4.])
Josef
<snip>
Chuck
_______________________________________________ SciPy-User mailing list SciPy-User@scipy.org http://mail.scipy.org/mailman/listinfo/scipy-user
On Sat, Feb 27, 2010 at 12:16 PM, <josef.pktd@gmail.com> wrote:
On Sat, Feb 27, 2010 at 2:07 PM, Charles R Harris <charlesr.harris@gmail.com> wrote:
On Sat, Feb 27, 2010 at 12:03 PM, Ivo Maljevic <ivo.maljevic@gmail.com> wrote:
Another suggestion, and this one also makes the parallel to Matlab. Make sure it works for vectors:
a=np.array([7,11,250]) nextpow2(a) array([ 3., 4., 8.])
In [27]: tab.searchsorted([7,11,250]) Out[27]: array([3, 4, 8])
In [28]: np.frexp([7,11,250])[1] Out[28]: array([3, 4, 8], dtype=int32)
It looks like frexp is doubling the array if n is already an integer power
np.frexp(2**np.arange(5)) (array([ 0.5, 0.5, 0.5, 0.5, 0.5]), array([1, 2, 3, 4, 5])) np.arange(5) array([0, 1, 2, 3, 4]) np.ceil(np.log2(2**np.arange(5))) array([ 0., 1., 2., 3., 4.])
Sure enough ;) Looks like the mantissa need to be checked for .5. In [22]: def nextpow2(x): m, e = np.frexp(np.abs(x)); return e - (m == .5) ....: In [23]: nextpow2(.5) Out[23]: -1 In [24]: nextpow2(.25) Out[24]: -2 In [25]: nextpow2(1) Out[25]: 0 In [26]: nextpow2(2) Out[26]: 1 Which brings up the question: is nextpow2(.5) 0 or -1? Chuck
-1 Which brings up the question: is nextpow2(.5) 0 or -1?
On Sat, Feb 27, 2010 at 4:14 PM, Ivo Maljevic <ivo.maljevic@gmail.com> wrote:
-1
Which brings up the question: is nextpow2(.5) 0 or -1?
frexp has a discontinuity at zero, same as matlab (or is this a computer science definition) matlab:
nextpow2(0.5) ans = -1 nextpow2(0.25) ans = -2
nextpow2(0.125) ans = -3 nextpow2(1e-300) ans = -996 nextpow2(0) ans = 0
Josef
_______________________________________________ SciPy-User mailing list SciPy-User@scipy.org http://mail.scipy.org/mailman/listinfo/scipy-user
On Sat, Feb 27, 2010 at 2:38 PM, <josef.pktd@gmail.com> wrote:
On Sat, Feb 27, 2010 at 4:14 PM, Ivo Maljevic <ivo.maljevic@gmail.com> wrote:
-1
Which brings up the question: is nextpow2(.5) 0 or -1?
frexp has a discontinuity at zero, same as matlab (or is this a computer science definition)
matlab:
nextpow2(0.5) ans = -1 nextpow2(0.25) ans = -2
nextpow2(0.125) ans = -3 nextpow2(1e-300) ans = -996 nextpow2(0) ans = 0
Yeah, zero is special, I was thinking -inf, but it's not an integer. Chuck
Thanks Robert. I figured that after I sent the question, but I am not sure if the search based solution offers any advantage in terms of speed. Besides, and not that it matters really, but that solution restricts the usage to max 2**32 integer result on 32 bit machines, whereas log based one doesn't:
def nextpow2(n): ... return 2**(np.ceil(np.log2(n))) ... np.iinfo(np.intp).bits 32 n=2**32+300.5 n 4294967596.5 nextpow2(n) 8589934592.0
On 27 February 2010 13:02, Robert Kern <robert.kern@gmail.com> wrote:
On Sat, Feb 27, 2010 at 04:02, Ivo Maljevic <ivo.maljevic@gmail.com> wrote:
David, Nice way of avoiding log2, but how do you determine the length of your tab array?
The input arrays can only be so large. Even on 64-bit machines, the table need not have more than 64 entries.
tab = 2 ** np.arange(np.iinfo(np.intp).bits)
-- Robert Kern
"I have come to believe that the whole world is an enigma, a harmless enigma that is made terrible by our own mad attempt to interpret it as though it had an underlying truth." -- Umberto Eco _______________________________________________ SciPy-User mailing list SciPy-User@scipy.org http://mail.scipy.org/mailman/listinfo/scipy-user
On Sat, Feb 27, 2010 at 12:32, Ivo Maljevic <ivo.maljevic@gmail.com> wrote:
Thanks Robert. I figured that after I sent the question, but I am not sure if the search based solution offers any advantage in terms of speed. Besides, and not that it matters really, but that solution restricts the usage to max 2**32 integer result on 32 bit machines, whereas log based one doesn't:
The point of nextpow2() is to find a good size array for doing the FFT. It is not a general purpose computational routine. On 32-bit machines, you cannot have an array larger than 2**32. -- Robert Kern "I have come to believe that the whole world is an enigma, a harmless enigma that is made terrible by our own mad attempt to interpret it as though it had an underlying truth." -- Umberto Eco
On Sat, Feb 27, 2010 at 11:32 AM, Ivo Maljevic <ivo.maljevic@gmail.com>wrote:
Thanks Robert. I figured that after I sent the question, but I am not sure if the search based solution offers any advantage in terms of speed. Besides, and not that it matters really, but that solution restricts the usage to max 2**32 integer result on 32 bit machines, whereas log based one doesn't:
def nextpow2(n): ... return 2**(np.ceil(np.log2(n))) ... np.iinfo(np.intp).bits 32 n=2**32+300.5 n 4294967596.5 nextpow2(n) 8589934592.0
For (ieee) floats In [17]: def nextpow2(x) : return np.ldexp(1., np.frexp(x)[1]) ....: In [18]: nextpow2(2**32+300.5) Out[18]: 8589934592.0 Chuck
participants (8)
-
Anne Archibald -
Charles R Harris -
David Cournapeau -
David Goldsmith -
Ivo Maljevic -
josef.pktd＠gmail.com -
Nils Wagner -
Robert Kern