Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Drastically speed up spectrums by zero-padding to a power-of-two. #6

Open
wants to merge 2 commits into
base: master
Choose a base branch
from

Conversation

kousu
Copy link

@kousu kousu commented Mar 21, 2016

Thank you for this code. It is very helpful to have all the methods in estimate_frequency.py laid out and contrasted.

I'm not totally sure about the change to thd.

This should be the numpy default, IMO, but what can you do?
@kousu
Copy link
Author

kousu commented Mar 21, 2016

..wait. This causes freq_from_hps() to give drastically wrong results. Why would that be? Is there a way to salvage it?

@endolith
Copy link
Owner

SciPy uses FFTPACK which is optimized for multiples of 2, 3, 5, which I made a function for here:

scipy/scipy#3144

The function is _next_regular https://github.com/endolith/scipy/blob/master/scipy/signal/signaltools.py#L246 but planning on changing it to _next_opt_len in the future: https://github.com/endolith/scipy/blob/czt/scipy/fftpack/helper.py#L49

@kousu
Copy link
Author

kousu commented Mar 22, 2016

Oh great. So this will this end up in numpy.fft.rfft() and then it will just always be fast? That would be excellent.

In the meantime, I don't suppose you have a clue why the padding I added is breaking your freq_from_hps()? I would like to use this function but it's too slow over a whole corpus. Using the longer fft means lower frequencies are captured, so I guess your warning about low frequency noise applies, but I don't know how to fix that. Should it be enough to manually clip the search range?

@kousu
Copy link
Author

kousu commented Mar 22, 2016

Oh, or are you saying that I should use _next_opt_len() instead of round2()?

@endolith
Copy link
Owner

So this will this end up in numpy.fft.rfft() and then it will just always be fast?

No, but czt can replace prime-length ffts in the future

Oh, or are you saying that I should use _next_opt_len() instead of round2()?

Yes, instead of powers of 2.

I don't suppose you have a clue why the padding I added is breaking your

What are you expecting to get and what are you actually getting?

…ound2()

Depending on the BLAS you have there are myriad optimal sizes. [For example](scipy/scipy#3144 (comment))
the optimal length under one implementation is any of 2^a 3^b 5^c 7^d 11^e 13^f,
where e+f is either 0 or 1, and the other exponents are arbitrary. It's just a
matter of how the implementation chooses to carve up subproblems. Halves---making
powers of 2 fast---is traditional, but not the only choice.
@kousu
Copy link
Author

kousu commented Mar 24, 2016

You know what? It was my sample. I was going to make a before and after test case and post it here, but as I waited and waited for the code to run across the vast majority of my corpus, I found that using a faster size or not only changes how the mistakes are made, not where. The mistakes almost always give different pitches, but the files that give mistakes trouble tend to be the same (there's 13 the unmolested version gets wrong where the round2()d version doesn't, and 8 vice versa). I think when I posted that I was tired and not looking that closely. I thought my Oboe corpus was drastically worse under rounding.

Thanks again for your code. I've updated the PR as you requested.

@endolith
Copy link
Owner

actually I'd be happier if next_regular were just copied into common.py as a public function instead of importing a private name that only exists in certain scipy versions and won't exist in future versions.

@endolith
Copy link
Owner

Now it's a public function next_fast_len.

@endolith
Copy link
Owner

I don't think this code will work as written. Have you tested it before and after with known frequencies? The new FFT length will be next_fast_len(N), but you're still using N to find the frequency in Hz: fs * i_interp / N. But i_interp will be shifted because the new spectrum is stretched out?

I am trying to clean up the project so I can't test this yet

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

2 participants