DSPRelated.com
Forums

Question about DSP in Matlab (beginner)

Started by JCO_DSP May 23, 2015
On 5/25/15 10:32 PM, Tim Wescott wrote:
> On Mon, 25 May 2015 07:53:20 -0500, Cedron wrote: >
...
>> >> Tim, I have to agree with Kaz on this one. When I first read your post, >> I thought it was unduly harsh and discouraging. Mathematics has often >> been described as syntactic encapsulation of patterns. Seeing the >> patterns first makes it way easier to understand the mathematical >> statement that describes them. That is in accordance with the >> Scientific Method, in which the first two steps are making observations >> and then formulating a hypothesis. In reverse, when you are learning >> what a mathematical statement means, it is generally quite useful to see >> an example of the pattern that the statement describes to understand it >> better. So to the OP, keep playing around, keep asking questions, and I >> hope you are enjoying learning. >> >> Speaking of comp.dsp negative nellies/bullies, Rune also appears to be >> gone.
i haven't heard from either, neither. Rune ... Vlad ... pussycats.
> > I'm sorry that you feel that way. Matlab will tell you that 1+1 = 2. > But it won't tell you that 1+1+1+1+... tends to infinity. If you don't > know the latter, then you don't know enough math to understand DSP theory.
to different degrees. there are a few "civilians" (often audio majors at some school) that might be able to understand some stuff by just tooling around MATLAB. but depth of understanding requires some formal math vocabulary. however, when i am at the dsp.stackexchange.com site, there are those who are far more mathematically "sophisticated" than me (they worry about ROC more than i do). i mean, it's been 35 years since grad school and i couldn't tell you shit about how to derive a Kalman filter anymore. bestest regrads (lest i be cast among the negative nellies/bullies of comp.dsp). -- r b-j rbj@audioimagination.com "Imagination is more important than knowledge."
On 5/25/15 10:32 PM, Cedron wrote:
> [...snip...] > > Vladimir V., is that you? Or just another angry dude?
no, no! Dmitry is not Vlad. but, if i recall, they have exchanged words. i think Vlad has picked on Dmitry.
> > Let me try to answer for the OP. > >> >> Why on Earth would anyone use Fourier transform to measure the frequency >> of a simple audio signal ?
and is this the frequency of a sinusoid or the fundamental frequency of a periodic (or quasi-periodic) audio signal? two sorta different problems.
> > Because it is a stepping stone to understanding how to deal with more > complicated signals. > >> And why do you think regular "experts" is this group are giving you >> correct answers ? >> > They are more likely to give consensus answers that are hopefully correct.
consensus won't always get you correctness (and often misses vision).
> If they don't, someone else will call them out. The resulting discourse > can be both entertaining and educational. > >> Google is your friend: >> http://www.google.com/patents/US7124075 >> http://www.gamedev.ru/flame/forum/?id4654 > > Patent abuse for obvious ideas are a whole another topic.
Dmitry's patent isn't totally obvious. i seem to remember the use of a histogram for, what is a similar output to autocorrelation (i said *similar*, Dmitry, not the same) to pick on pitch candidates. but i still think it's a sorta jazzed-up AMDF with a few non-linear functions in between, one of them was a step function. and i am dubious about any step functions because information is lost and that means you can construct a signal that will fool the pitch detector and you could get octave errors or the like. -- r b-j rbj@audioimagination.com "Imagination is more important than knowledge."
On 5/25/15 10:37 PM, Tim Wescott wrote:
>
...
> Aside from the absurd cost model, Matlab
i have, just very recently, discovered a "cheap" MATLAB that wasn't a student version. won't say who told me. http://www.mathworks.com/pricing-licensing/?intendeduse=home seems to not be handicapped in any manner.
> is a fine tool for experimenting > to make sure that your theoretical constructs do what you expect.
i might use it for other purposes, in addition.
> I use > Scilab for exactly that -- which is your "Matlab as a learning tool". > > Matlab (or Scilab) is also a fine tool for crunching lots of numbers > _after_ you have used theoretically sound methods to determine the > details of the number crunching algorithms. > > But if you think that typing > > y = fft(x); > plot(abs(y)) > > means that you thoroughly understand the FFT, then you've been mislead.
that's fer sherr. -- r b-j rbj@audioimagination.com "Imagination is more important than knowledge."
>On Mon, 25 May 2015 07:53:20 -0500, Cedron wrote: > >>>>On Sat, 23 May 2015 12:10:22 -0500, JCO_DSP wrote: >>>> >>>>> Hi, >>>>> >>>>> I'm starting to learn DSP (Matlab ) >>>> >>>>If you are equating knowing Matlab with knowing DSP, then your error >>>is >>>> >>>>insurmountable until you understand that DSP is a thing you do, >>>>mathematics is the tool you use to understand it, and Matlab is a >tool >>> >>>>you use, _after_ you understand the math, to do the math. >>>> >>>>-- >>>>www.wescottdesign.com >>>
[...snip...] My this is an active thread. My disagreement with the above statement is the implication of precedence order. As I've elaborated in other responses, playing around with the numbers to see how they work can give a familiarity that will make the math easier to learn. I do not think it is necessary to know the math first. On the other hand, it is not necessary to play around first either.
>I'm sorry that you feel that way. Matlab will tell you that 1+1 = 2. >But it won't tell you that 1+1+1+1+... tends to infinity. If you don't > >know the latter, then you don't know enough math to understand DSP >theory. > >-- >www.wescottdesign.com
Well, I've never used Matlab, but I suspect it wouldn't do very well with "come up with a function that maps from R to RxR that is completely space filling" either. I learned that 1/3 + 1/9 + 1/27 + 1/81 approached 1/2 by playing on a calculator way before I learned about infinite series in a formal way. I'm not sure what your are feeling sorry about which way I feel. If it was the bully remark, that was referring to Vlad and Rune, not you. I still think the timbre of your remark was discouraging, but far from being bullying. Now, if we were talking about using a chain saw, a little formal training before playing with it might be a really good idea. Ced --------------------------------------- Posted through http://www.DSPRelated.com
robert bristow-johnson  <rbj@audioimagination.com> wrote:

>Dmitry's patent isn't totally obvious. [...]
I can tell you firsthand that you can invent a novel, non-obvious invention that is actually original and creative, but by the time the patent lawyers get done with filing it, it will look mostly completely obvious, or at least 95% of "what is claimed" is totally obvious. Somewhere down in there may be an invention; but you'll have to dig deep. Steve
On 5/25/15 3:19 PM, Cedron wrote:
...
>> Repeating: >>> Unless you can derive the equations in my blog article titled "DFT Bin >>> Value Formulas for Pure Real Tones" (at dsprelated dot com), you don't >>> understand leakage mathematically. Sorry, you don't. >> >> I think that's an incredibly arrogant statement. >> >> And, btw, other people have done derivations to the same purpose as >> your article and come up with cleaner results. >> > > Thank you. I enjoy writing in the "authoritative voice". > > If you don't mind, could you provide some links for these derivations? > > rbj thought that the exact frequency formula had been solved, yet I showed > him links to your frequency estimation page and the dspguru page, there > was no such formula to be found, just estimators. No reply.
okay, i only barely remember any of this. maybe i blew it off, Ced. sorry. first of all, regarding the correspondence i had with Eric that is alluded to at http://www.ericjacobsen.org/fe.htm , i remember there was a factor of 2 that he had that i did not (or the other way around) that i never grok'd what was the source of our discrepency. also, the quadratic thing that i was referring to was not a particularly sophisticated thing and *is* an estimate, with one exception that i am aware of: if you use a gaussian window and are looking at log magnitude of the FFT, then for a single sinusoidal component it *should* fit the quadratic exactly (if leakage from the negative frequency component is not a problem). now, looking at your blog Ced, (first lemme say i don't like your notation, but i am a firm believer in the notational convention of Scripture, the Word of O&S, so you got your k's and n's mixed up and that requires effort from me to keep straight). now i know what alpha and beta_n are. can i take it that your formula solves for alpha in a closed form from Eq. 17? alpha appears 4 places in that equation. do you have an expression for alpha in terms of all of the Z_n (which are the bins of the FFT)? what is it? and, seems to me, that this analysis is dependent on a rectangular window (or, as Eric and i debated 20 years ago, no window at all). i wanna see a nice clear expression of alpha in terms of the Z_n bins in order for me to understand how you solved the problem, Ced. BTW, the way i think it's a "solved problem" is in the sense of this old paper https://secure.aes.org/forum/pubs/journal/?elib=7044 and i don't remember how it was all done, but at one time i think i did know. it wasn't cheap and i think the authors had to try a bunch of different frequencies and see which fit best. or maybe what i was thinking of was "estimating" (if that's the word you wanna use) the frequency, phase, frequency sweep-rate, amplitude, and amplitude ramp-rate of a frame of audio that was windowed off with a gaussian window. *if* you assume no leakage from other components (or the negative frequency component) because of the rapid attenuating skirt of the gaussian window in the frequency domain, then you can "estimate" those parameters pretty dang directly and i wrote a paper about it in 2001: https://www.researchgate.net/publication/3927319_Intraframe_time-scaling_of_nonstationary_sinusoids_within_the_phase_vocoder . don't squeal to IEEE that i put it up on ResearchGate. it's the only IEEE pub i had ever done (and don't expect to do another one so i don't give a rat's ass what IEEE will do to me). L8r, -- r b-j rbj@audioimagination.com "Imagination is more important than knowledge."
On Monday, May 25, 2015 at 11:59:17 PM UTC-4, robert bristow-johnson wrote:
> On 5/25/15 3:19 PM, Cedron wrote: > ... > >> Repeating: > >>> Unless you can derive the equations in my blog article titled "DFT Bin > >>> Value Formulas for Pure Real Tones" (at dsprelated dot com), you don't > >>> understand leakage mathematically. Sorry, you don't. > >> > >> I think that's an incredibly arrogant statement. > >> > >> And, btw, other people have done derivations to the same purpose as > >> your article and come up with cleaner results. > >> > > > > Thank you. I enjoy writing in the "authoritative voice". > > > > If you don't mind, could you provide some links for these derivations? > > > > rbj thought that the exact frequency formula had been solved, yet I showed > > him links to your frequency estimation page and the dspguru page, there > > was no such formula to be found, just estimators. No reply. > > > okay, i only barely remember any of this. maybe i blew it off, Ced. sorry. > > first of all, regarding the correspondence i had with Eric that is > alluded to at http://www.ericjacobsen.org/fe.htm , i remember there was > a factor of 2 that he had that i did not (or the other way around) that > i never grok'd what was the source of our discrepency. also, the > quadratic thing that i was referring to was not a particularly > sophisticated thing and *is* an estimate, with one exception that i am > aware of: if you use a gaussian window and are looking at log magnitude > of the FFT, then for a single sinusoidal component it *should* fit the > quadratic exactly (if leakage from the negative frequency component is > not a problem). > > now, looking at your blog Ced, (first lemme say i don't like your > notation, but i am a firm believer in the notational convention of > Scripture, the Word of O&S, so you got your k's and n's mixed up and > that requires effort from me to keep straight). now i know what alpha > and beta_n are. can i take it that your formula solves for alpha in a > closed form from Eq. 17? alpha appears 4 places in that equation. do > you have an expression for alpha in terms of all of the Z_n (which are > the bins of the FFT)? what is it? and, seems to me, that this analysis > is dependent on a rectangular window (or, as Eric and i debated 20 years > ago, no window at all). > > i wanna see a nice clear expression of alpha in terms of the Z_n bins in > order for me to understand how you solved the problem, Ced. > > BTW, the way i think it's a "solved problem" is in the sense of this old > paper https://secure.aes.org/forum/pubs/journal/?elib=7044 and i don't > remember how it was all done, but at one time i think i did know. it > wasn't cheap and i think the authors had to try a bunch of different > frequencies and see which fit best. > > or maybe what i was thinking of was "estimating" (if that's the word you > wanna use) the frequency, phase, frequency sweep-rate, amplitude, and > amplitude ramp-rate of a frame of audio that was windowed off with a > gaussian window. *if* you assume no leakage from other components (or > the negative frequency component) because of the rapid attenuating skirt > of the gaussian window in the frequency domain, then you can "estimate" > those parameters pretty dang directly and i wrote a paper about it in > 2001: > https://www.researchgate.net/publication/3927319_Intraframe_time-scaling_of_nonstationary_sinusoids_within_the_phase_vocoder > . don't squeal to IEEE that i put it up on ResearchGate. it's the only > IEEE pub i had ever done (and don't expect to do another one so i don't > give a rat's ass what IEEE will do to me). > > L8r, > > > -- > > r b-j rbj@audioimagination.com > > "Imagination is more important than knowledge."
Just to add to what robert wrote: Frequency estimation of a sinusoid is a very well studied field, and there is a huge amount of literature available, such as Gaussian interpolation and Gaussian ratio: https://ccrma.stanford.edu/~jos/sasp/Quadratic_Interpolation_Spectral_Peaks.html (as noted in the above: "the Gaussian window transform magnitude is precisely a parabola on a dB scale. As a result, quadratic spectral peak interpolation is exact under the Gaussian window." - the authors' book also has other non-Gaussian estimation algorithms) http://www.edn.com/design/analog/4352503/EDN--03-03-94-Ratio-detection-precisely-characterizes-signals-amplitude-and-frequenc (the link has a few problems, but the gist is that you can obtain very precise frequency estimates - I've used this technique a lot since I found out about it in the early 1990's - 12 or more decimal digits of accuracy is not at all unusual). http://www.ptodirect.com/Results/Patents?query=PN/5214708 (a patent related to the EDN article - and like the EDN article, it too has a few problems) In addition, a related Gaussian technique is used at CERN (home of the $10 billion LHC). The authors of the technique have published a couple of scientific papers. On of them also has a PhD dealing with spectral peak estimation, which is available on-line (for those folks who think that an error bound of 10 to the minus 12 is waaaay too imprecise). In addition, there are many other techniques. But I'm not going to search for the links right now. Maybe tomorrow. Kevin McGee
[...snip...]
>> >> rbj thought that the exact frequency formula had been solved, yet I >showed >> him links to your frequency estimation page and the dspguru page, >there >> was no such formula to be found, just estimators. No reply. > > >okay, i only barely remember any of this. maybe i blew it off, Ced. >sorry. >
No worries. If you had asserted that there was, I would have pressed for an answer. Saying that you thought there was is very non-committal, and in my mind, doesn't require justification. Unlike EJ saying there are cleaner versions of my bin equations and then not providing a cite.
>first of all, regarding the correspondence i had with Eric that is >alluded to at http://www.ericjacobsen.org/fe.htm , i remember there was > >a factor of 2 that he had that i did not (or the other way around) that > >i never grok'd what was the source of our discrepency. also, the >quadratic thing that i was referring to was not a particularly >sophisticated thing and *is* an estimate, with one exception that i am >aware of: if you use a gaussian window and are looking at log magnitude > >of the FFT, then for a single sinusoidal component it *should* fit the >quadratic exactly (if leakage from the negative frequency component is >not a problem). >
I was looking mainly at fe2.htm.
>now, looking at your blog Ced, (first lemme say i don't like your >notation, but i am a firm believer in the notational convention of >Scripture, the Word of O&S, so you got your k's and n's mixed up and >that requires effort from me to keep straight).
Sorry, those are the conventions I have been using for a long time. I developed my own stuff from scratch. I don't know who O&S is. Do you think it is worth editing all my blogs to change this?
>now i know what alpha >and beta_n are. can i take it that your formula solves for alpha in a >closed form from Eq. 17? alpha appears 4 places in that equation.
No, alpha only appears once in (17), on the left hand side. It appears twice in (16). I had to put a line break in some of the equations because of the limited display width.
>do >you have an expression for alpha in terms of all of the Z_n (which are >the bins of the FFT)? what is it? and, seems to me, that this analysis > >is dependent on a rectangular window (or, as Eric and i debated 20 years > >ago, no window at all).
Yes. I prefer the "no window at all" point of view. As I explain later, the nature of the weighting vector makes it as though there is a virtual near VonHann window being applied.
> >i wanna see a nice clear expression of alpha in terms of the Z_n bins in
That would be (17), (18), and (19). Only three bins are used. I solved the equations for three adjacent bins, but any three or more bins could be used. Also, the formula is not dependent on having the bin set located at the peak. In the numerical examples that follow, I calculate the frequency from the peak, the DC bin, and the Nyquist bin. R1 = exp( -i 2 pi / N ) ~=~ 1
> >order for me to understand how you solved the problem, Ced. >
That's what the derivation is supposed to be for. Conceptually, I solved the inverse of the bin value formula. It is "clean" enough to be invertible in closed form. The bin values are a function of the frequency (and phase/magnitude) so the frequency is then a function of the bin values.
>BTW, the way i think it's a "solved problem" is in the sense of this old > >paper https://secure.aes.org/forum/pubs/journal/?elibp44 and i don't >remember how it was all done, but at one time i think i did know. it >wasn't cheap and i think the authors had to try a bunch of different >frequencies and see which fit best. > >or maybe what i was thinking of was "estimating" (if that's the word you >
A numerical solution is not the same as a closed form solution. The abstract alludes to an iterative approach. The paper also seems to be dealing with a varying signal, my formula pertains to a pure unwavering tone, so its significance is more theoretical than practical.
>wanna use) the frequency, phase, frequency sweep-rate, amplitude, and >amplitude ramp-rate of a frame of audio that was windowed off with a >gaussian window. *if* you assume no leakage from other components (or >the negative frequency component) because of the rapid attenuating skirt > >of the gaussian window in the frequency domain, then you can "estimate" > >those parameters pretty dang directly and i wrote a paper about it in >2001: >https://www.researchgate.net/publication/3927319_Intraframe_time-scaling_of_nonstationary_sinusoids_within_the_phase_vocoder > >. don't squeal to IEEE that i put it up on ResearchGate. it's the only >
I won't. But ResearchGate wants me to register to see it and I don't want to do that.
>IEEE pub i had ever done (and don't expect to do another one so i don't > >give a rat's ass what IEEE will do to me). > >L8r, > > >-- > >r b-j rbj@audioimagination.com > >"Imagination is more important than knowledge."
The DFT is lossless. That is why I think the terms "leakage" and "uncertainty principle" are so inappropriate. Closely spaced tones can be resolved in the time domain by plotting their envelope (if they are the only two tones in the signal). They can also be resolved from the DFT by using a follow the gradient type iterative method using the bin value equations to produce basis vectors, sort of similar to my approach of phase and magnitude determination (not estimation) that I document in my latest blog. Thanks for your reply. I'd be happy to answer any questions on any of the derivations you may have. I would also recommend that you check out my first two articles. Ced --------------------------------------- Posted through http://www.DSPRelated.com
kevinjmcee  <kevinjmcgee@netscape.net> wrote:

>In addition, a related Gaussian technique is used at CERN (home of the >$10 billion LHC). The authors of the technique have published a couple >of scientific papers. On of them also has a PhD dealing with spectral >peak estimation, which is available on-line (for those folks who think >that an error bound of 10 to the minus 12 is waaaay too imprecise).
Q. What's the difference between a Physicist and an Engineer? A. Physicist: "If we apply this cool Mexican Hat filter to our dataset, we can identify spectral features that haven't been seen since the Big Bang!!" Engineer: "Don't tell management, but everybody knows that when you want to produce spectral features that aren't actually present in your data, try a Mexican Hat filter." S.
On Mon, 25 May 2015 22:51:34 -0400, robert bristow-johnson wrote:

> On 5/25/15 10:32 PM, Tim Wescott wrote: >> On Mon, 25 May 2015 07:53:20 -0500, Cedron wrote: >> > ... >>> >>> Tim, I have to agree with Kaz on this one. When I first read your >>> post, >>> I thought it was unduly harsh and discouraging. Mathematics has often >>> been described as syntactic encapsulation of patterns. Seeing the >>> patterns first makes it way easier to understand the mathematical >>> statement that describes them. That is in accordance with the >>> Scientific Method, in which the first two steps are making >>> observations and then formulating a hypothesis. In reverse, when you >>> are learning what a mathematical statement means, it is generally >>> quite useful to see an example of the pattern that the statement >>> describes to understand it better. So to the OP, keep playing around, >>> keep asking questions, and I hope you are enjoying learning. >>> >>> Speaking of comp.dsp negative nellies/bullies, Rune also appears to be >>> gone. > > i haven't heard from either, neither. Rune ... Vlad ... pussycats. > > >> I'm sorry that you feel that way. Matlab will tell you that 1+1 = 2. >> But it won't tell you that 1+1+1+1+... tends to infinity. If you don't >> know the latter, then you don't know enough math to understand DSP >> theory. > > to different degrees. there are a few "civilians" (often audio majors > at some school) that might be able to understand some stuff by just > tooling around MATLAB. but depth of understanding requires some formal > math vocabulary. however, when i am at the dsp.stackexchange.com site, > there are those who are far more mathematically "sophisticated" than me > (they worry about ROC more than i do). i mean, it's been 35 years since > grad school and i couldn't tell you shit about how to derive a Kalman > filter anymore. > > bestest regrads (lest i be cast among the negative nellies/bullies of > comp.dsp).
Oh, deriving a Kalman filter is easy -- just pick up a good Kalman filter book, turn to the appropriate chapter (usually 2, although sometimes it'll be in 1 or 3), and start reading... -- www.wescottdesign.com