Wouldn't it be more efficient and accurate to just sum up the product of the voltage times the current? What am I missing?
Perhaps, this whole power factor thing is a bit of a mystery to me.
But I'm thinking:
1) If the load is purely resistive then the voltage and current are in phase and the average of the product of the voltage times the current will give some positive value.
2) If the load is inductive the current will lag behind the driving voltage. Then the average of the product of the voltage times the current will be some lesser value.
3) If the load is capacitive the current will lead the voltage. Then the average of the product of the voltage times the current will be some lesser value.
4) In the extreme the phase shift is 90 degrees. Then the resulting average is zero.
I guess if you have measured the current and voltage levels as you go, then those lesser values in 1) and 2) give you an idea of the phase shift and hence power factor some how.
I was thinking it would be an idea just to:
1) Take the average of the voltage multiplied by a reference sine wave at the required frequency, Vs.
2) Take the average of the voltage multiplied by a reference cosine wave at the required frequency, Vc
3) Do the same for current giving Is and Ic.
Now arctan(Vs / Vc) gives the phase of the incoming voltage relative to the internal refence signals.
Similarly arctan(Is / Ic) gives the phase of the incoming current signal.
The difference between the two is the phase shift of current vs voltage.
Your technique will work fine as long as the voltage and current are pure sine waves. There is a problem if there are harmonics, as you mentioned earlier. You can put the current in the imaginary part, but it must be scaled the same as the voltage. Also, I think the resulting phase is half of the actual phase difference, but I might be wrong on that.
Intantaneous power is defined as P = I * V. This needs to be integrated (or summed) over a complete cycle to get the average value for AC power. The power factor is the ratio of the actual power versus the power if I and V are in phase. I think the in-phase power would be the product of the RMS values of V and I.
For now I'm just paying attention to the primary frequency, 60 Hz. Later I may try working with harmonics, but since this is for an energy meter for a whole house, it probably isn't necessary.
Here's the little bit of code I had to write using your awesome FFT object:
Thanks for the code. One question I had, does it assume a sampling rate of 1024 Hz, or is it just a coincidence that the result locations line up with their frequencies? I've done plenty of FFT in Sysquake and MATLAB, and have never noticed this before.
So what you are saying is that given a bunch of samples of current and voltage, 1024 per cycle each say, we can get the power factor by:
1) Calculating the RMS values of the voltage and current samples, Vrms and Irms.
2) Get RMS power from that: Prms = Vrms * Irms
3) Now calculate the instantaneous power at each sample point Vn * In
4) Take the average of that to get the real power Pr
The Power Factor = Pr / Prms
Further that this is correct even if the waves forms are distorted and full of harmonics.
.
This rather makes the FFT a big hammer to crack a nut:)
...does it [the FFT] assume a sampling rate of 1024 Hz, or is it just a coincidence that the result locations line up with their frequencies?
Well if you look at it you see that the FFT has no idea about time and no idea about your Hz.
When you see an amplitude in an output bin, bx[60] for example, it just means "60 cycles per 1024 samples". If you happen to be sampling at 1024 samples per second then that represents 60Hz in reality. If you were sampling at 16384 samples per second then bx[60] would represent 960Hz. And so on. I think I have seen this referred to as the "Fourier Frequency" in some literature. It's the frequency within the sample set rather than any real frequency.
For now I'm just paying attention to the primary frequency, 60 Hz. Later I may try working with harmonics, but since this is for an energy meter for a whole house, it probably isn't necessary.
Turns out that FFT might be overkill for measuring power factor. But seeing as you have it working nicely I think you should consider calculating the percentage harmonic distortion from the data you have available in the FFT output.
We might guess that the voltage distortion is quite low and constant but the current distortion could be large if there are non-linear loads like rectifiers or SCR controlled light dimmers.
So if you see a high distortion you know it's time to be looking at those harmonics for power factor or consider a different technique.
Actually I would have found that useful to have a few years back when we had trouble getting some equipment working on an UPS. Turned out it was monitoring mains frequency and amplitude for safety reasons but the UPS was generating a lot of harmonics that confused it causing it to trip out all the time.
Anyway, I have now learned more about the mysterious power factor than I ever though I would:)
we can get the power factor by:
1) Calculating the RMS values of the voltage and current samples, Vrms and Irms.
2) Get RMS power from that: Prms = Vrms * Irms
3) Now calculate the instantaneous power at each sample point Vn * In
4) Take the average of that to get the real power Pr
The Power Factor = Pr / Prms
Further that this is correct even if the waves forms are distorted and full of harmonics.
Heater,
Your description is much better than the way I was trying to explain it. And, this works with any type of waveform -- not only for pure sinusoids.
Heater, another thing I forget to mention. My voltage and current sense circuits are referenced to a virtual ground, for input to an ADC. Using your code, I can easily get the offsets by using the magnitude of the 0 index FFT result. Maybe it is a big hammer to crack a nut, but it works and is making all parts of my project easier. Since this is for an energy meter, and I read EU charges for certain THD levels, it'll help a lot for that also.
Jay,
What do you use for a current sense transformer?
RS_Jim
CR8420-1000-G from CR Magnetics. But any should work, as long as you adjust the burden resistance too keep the voltage drop linear with primary current. I'll be posting all the details of the project very soon.
I'm a student using the Propeller for an EE group design project in school. The group is building a 2 octave keyboard synthesizer. Here is a picture of the keyboard/keyboard frame:
. Now we're working on the synthesis algorithms. The final project needs to be able to produce square, sawtooth, and triangle waveforms (based on a switch on the board) Variable pressure sensors have been implemented to detect the pressure being applied to each key, and 4 MCP3208's are used to read in the voltage value that those sensors produce into the Propeller IC. The FFT and the inverse FFT are essentially the same, except the inverse fft samples produced are divideded by the sample size N, so we would ideally like to directly synthesize the frequencies needed based on the frequency of the specific key (keys in the future) being pressed. One problem I have found is that, with a 1024 point algorithm, the resolution is not fine enough to get the exact frequencies needed for each key (since our required bandwidth is so high because of the harmonics of the triangle, sawtooth, and triangle waveforms). I know the resolution can't be fine enough to hit all the frequencies for synthesis perfectly, but I'd like to modify the Heater FFT algorithm to at least a 2048 point algorithm, and possibly higher, so that the frequencies synthesized could be closer to the actual ideal frequencies of the keys. The FFT file says that higher point FFTs may require modification of the 'twiddle factor tables' and the 'fixed point format'. I was hoping to get an explanation for these modifications so that I could attempt to do them. Any help or advice regarding this (or the project in general) would be really appreciated. Thanks in advance!
It seems to me to be extreme Overkill; using an FFT to produce these three simple waveforms at the frequencies needed. Would it not be much more practical to use much simpler to go high to low with a counter for the square wave, and use an incrementing/decrementing counter at some multiple of the freq to produce a triangle wave, then modify that to have a sharp leading edge with a decrementing counter for the saw tooth wave? I know that is how I do produce these waveforms. Does this project require that you use a FT?
It may be....but since our project does require that we can produce the signals in an envelope, and we would like to be able to press multiple keys at a time, it seems to me like directly synthesizing the frequencies from the frequency domain would be the most practical way of getting the job done. We could implement the signal envelope by changing the magnitude of the frequencies being generated from the frequency domain in real time.
To produce your envelope just mod your amplitude by another series of directly generated values. As for mixing, if you output through a resister DAC at 512 KHz to produce the wave, you can simply interleave up to 8 channels that once filtered with an appropriate capacitor would produce a very smooth mixing for the waveform. Do they also teach you to use a 20 ton hydraulic press to put a thumb tack in cork board?
Take your unmodified output value and divide it by the number of steps from 0 output to max output (i recommend using the high bits to drive the DAC to simplify this) (this should be a power of two making for a simple shift), then multiply it by the corresponding value in the generated envelope wave. It does not get much simpler than that.
This is more brief than originally intended. Firefox crashed on me 4 times while I was typing.
Wow I really appreciate the insight. We were obviously headed in the wrong direction. I just have a few more quick questions and this will be the last post on this thread because I know we're getting off topic. You said:
As for mixing, if you output through a resister DAC at 512 KHz to produce the wave, you can simply interleave up to 8 channels
When you say 'interleave', what do you mean by that? Do you mean add the modulated signals together before we output them to the resistor DAC? And I'm sure I'm going to get slammed for this question but when you say:
Take your unmodified output value and divide it by the number of steps from 0 output to max output (i recommend using the high bits to drive the DAC to simplify this)
Does this mean that the envelope range should go from 0 to 2^n, so that when the modded signal it sent to the dac it effectively ranges from 0 to 1? Thanks a lot for all your help, this probably saved our project.
When you say 'interleave', what do you mean by that? Do you mean add the modulated signals together before we output them to the resistor DAC?
No, I mean put one channels current on the output pins, fallowed by the next, repeat in a round robin fasion for each active channel. becouse this is done at quite a bit higher freq than the range of audibly audio, once you filter the final analog through a ~ 0.1 micro-farad cap you end up with a very smooth mixer, no math needed. You can even scope the output and you could not tell that it was doen with this method, it will appear properly mixed.
Does this mean that the envelope range should go from 0 to 2^n, so that when the modded signal it sent to the dac it effectively ranges from 0 to 1? Thanks a lot for all your help, this probably saved our project.
You could do that. Assuming you are using an 8-bit resistor dac (for the purpose of example), I would have the resister DAC correspond to bits 24-31 of your working value, and thus the range would be from 0 to $FF000000 hex in steps of $01000000, this will allow the math to be done with reasonable ease WITH OUT needing a fixed point or floating point lib.
Not sure if this FFT is what you need but a few comments anyway,
The FFT and the inverse FFT are essentially the same, except the inverse FFT samples produced are divideded by the sample size N"
Depends. According to the definition on wikipedia the FFT has no 1/N but the inverse does. Anyway the heater_fft butterfly code does not have 1/N that is left until the magnitude calculation stage.
More importantly the inverse FFT has a different sign for the twiddle factors.
..so we would ideally like to directly synthesize the frequencies needed based on the frequency of the specific key
Not a bad idea. I fear the problem here is one of speed. What frequency are you hoping to go up to?
One problem I have found is that, with a 1024 point algorithm, the resolution is not fine enough to get the exact frequencies needed for each key
That is a problem. If you want a 1Hz resolution up to 20KHz you are going to need 40K samples per second. What frequency resolution do you think you need?
The FFT file says that higher point FFTs may require modification of the 'twiddle factor tables' and the 'fixed point format'. I was hoping to get an explanation for these modifications so that I
could attempt to do them
The "twiddle factor" tables were simply one half cycle of cos and minus sine (wx and wy respectively). They are scaled to +/- 4095. So to double the FFT size we double the size of those tables and recalculate those.
BUT the latest versions of heater_FFT don't define their own sin/cos tables they look up what the need in the Props PROM tables.
My worry is, and I have not worked this through much, that when you double the size of the FFT you also double the max values in the results. Things might start to overflow.
As I said the twiddle factor tables have a range +/-4095. You will see some shifts by 12 in the code where the multiplies are taking care of this fixed point scaling. On could reduce the twiddle
factor range to +/- 2047 or such and reduce the shifts accordingly to reduce overflow problems.
Now if you decide to generate your waves by a simpler route you might like to think about Direct Digital Synthysis(DDS). Here is a chap generating nice sine waves on an AVR with only half a dozen lines of AVR assembler code: http://www.myplace.nu/avr/minidds/index.htm
No, I mean put one channels current on the output pins, fallowed by the next, repeat in a round robin fasion for each active channel. becouse this is done at quite a bit higher freq than the range of audibly audio, once you filter the final analog through a ~ 0.1 micro-farad cap you end up with a very smooth mixer, no math needed. You can even scope the output and you could not tell that it was doen with this method, it will appear properly mixed.
Wow that sounds like a great solution! We will certainly check it out in the oscilloscope, and check the frequencies that it reads. I'm going to work on this and let everyone know how it goes. Thanks for all the help
Of course while you may use simpler methods than the FFT to generate the tones you want you could always add the FFT to the system and display a graphic of the tones being generated as you play.
Heater:
Good thought. The reason I had suggested the methods mentioned above is that using an reverse FT would take enough processing time that the effective sample rate would go through the floor, and while mixing could be accomplished, with only 20MIPs per core I have found that the quality is lost. With the Simpler methods you get just as good of mixing as you would using a reverse FT on a system fast enough to be capable of keeping ahead of the sample rate, with out any loss. and you still have enough time left over to communicate between COGs, thus the sky is the limit. Though I can definitely see showing the wave form on a display, and a reverse FT would be one way to accomplish this as the time constraints are not as important.
I agree. It's doubtful that the FFT would have the performance required to do the actual synthesis. Unless of course one splits the FFT over 2 or 4 COGs. But then it's using all your resources and probably quite a pain to get working nicely. Things get worse if you want to double or quadrupal the FFT size. which will more than double or qudrupal it's execution time.
My suggestion for using the FFT to display the spectrum of what is being synthesized is just icing on the cake for this project. As you say the time constraints are not there. Just sample, transform and display at whatever speed it goes around at.
P.S. I do like the style of that keyboard in the picture.
added a command that applies the Hamming window to input
added an option to calculate power instead of magnitude
With this new code, if you need to calculate magnitude, pass CMD_POWER | CMD_MAGNITUDE. This works exactly same as specifying CMD_MAGNITUDE in original code. If you only pass CMD_POWER, then no square root operation is performed, and on exit bx contains power of each freq component (power is magnitute squared).
Specifying CMD_HAMMING applies the Hamming window to input. If CMD_COMPLEX is also specified, both bx and by are affected by window function, otherwise only bx is affected. Hamming window only works for FFT_SIZE of 8192 and below.
Brilliant, this is a great contribution. I'll roll it into a version 3.0 of heater_fft after I have played with it for a bit.
In case anyone is wondering what this "Hamming Window" is about I think I can summarize it like this:
The FFT works very well if all the frequencies present in the input signal are such that their wavelengths fit exactly into the sample buffer, that is one cycle of a sine wave, say, starting at zero in the first sample, rising to a max, passing through zero in the middle of the buffer, going negative in the second half of the buffer and then returning to zero at the last sample position. Or two, three, four, etc cycles up to the nyquist frequency. With these "neat" waveforms in the input there are no discontinuities in amplitude or slope at the ends of the buffer and a nice clean spectrum is produced by the FFT.
BUT what happens with real world signals where you may have 10.3 cycles in the input? The waves don't fit neatly into the buffer, they don't "wrap around" nicely. There are discontinuities at the ends. These discontinuities appear like noise to the FFT which will now produce an output showing noise frequencies through out the spectrum.
The way to reduce that noise is to reduce the effects of the samples at the ends of the buffers, where the discontinuities occur. This is done by multiplying all the input samples by a function that is zero (or near zero) at the ends of the buffer but rises to a non-zero value at the center of the buffer. There are many such functions. A simple triangular function would work for example. Hamming is such a function based on half a cycle of cosine.
See wikkipedia for more info on this technique.
Now, I have not had much time to check it out yet but I did apply the Hamming window to the test tones that are set up in the heater_fft demo. Thing is with those tones I found when plotting the results that there was not much difference in output. The major difference being that all components are reduced to about 30% amplitude compared to not using the window.
I'm guessing things would look different when looking at lower frequencies and especially fractional frequencies that don't fit as whole cycles into the input buffers.
Andrey, do you have an application in mind for this, I'm sure we'd all like to hear about it.
P.S. This code look really nice in PZST:)
P.P.S. I'm a bit worried about this idea of growing the FFT size passed 1024 samples. There is, I think, a problem that as the FFT size grows the calculations are going to overflow sooner.
The major difference being that all components are reduced to about 30% amplitude compared to not using the window.
No surprise here. Each window function is characterized with coherent power gain. It measures the reduction in signal power due to the window function suppressing a coherent signal at the ends of the measurement interval. For Hamming window, coherent power gain is 0.54. Since you are measuring magnitude, not power, expect a loss of 1 - sqrt(0.54) ≈ 0,265
Andrey, do you have an application in mind for this, I'm sure we'd all like to hear about it.
I am just playing with various DSP techniques. Using your FFT, I think I could build an audio spectrum analyzer with unique features - 2 independent channels, 31 Hz resolution, virtually any number of frequency channels. And, most exciting - same unit can be used with virtually any display size. Add few shift registers - and make display on 8x8 LED matrix, or 32x16 matrix. Or hook up LCD panel (or even two of them ) for a giant DJ spectrum display.
Well I'm not surprised at the reduction in output. After all we are chopping down the useful signal at the ends as well as the disturbing discontinuities.
When you say "coherent signal" do you mean a signal that "fits" nicely in the sample buffer? Sounds like it.
I have also had visions of building such an audio spectrum analyzer or DJ "light organ" from FFT. I just don't have time to keep up with the endless stream of projects that present themselves:) It would be great to see.
There are those who suggest the Hartley Transform is the way to go. Supposedly faster and using less memory. As I don't begin to understand the thing I have not looked at Hartley codes much. It's taken me decades to finally realize how the FFT works:) There is a Hartley Transform for the Prop somewhere.
I'm going to try and create some non-coherent test input for heater_fft that shows off what the Hamming Window can do. Hopefully I can post some plots here soon.
Be aware that the execution time of heater_fft is dependent on input data. That's because the multiply routine has a speed dependent on the size of the numbers being multiplied. Perhaps a fixed time multiply loop would be better.
I'm not sure what I'm looking at there but it does not immediately look right.
I would have thought that a sine wave multiplied by the Hamming function would be alternating equally positive and negative, as it normally does, but with higher amplitude in the middle of the buffer. Your second plot shows a sine wave with a huge DC bias, no negative going.
Yes, it has DC offset large enough. This was done intentionally, so the resulting waveform envelope looks like the window itself. Here is how it looks for a waveform without DC offset:
Comments
Perhaps, this whole power factor thing is a bit of a mystery to me.
But I'm thinking:
1) If the load is purely resistive then the voltage and current are in phase and the average of the product of the voltage times the current will give some positive value.
2) If the load is inductive the current will lag behind the driving voltage. Then the average of the product of the voltage times the current will be some lesser value.
3) If the load is capacitive the current will lead the voltage. Then the average of the product of the voltage times the current will be some lesser value.
4) In the extreme the phase shift is 90 degrees. Then the resulting average is zero.
I guess if you have measured the current and voltage levels as you go, then those lesser values in 1) and 2) give you an idea of the phase shift and hence power factor some how.
I was thinking it would be an idea just to:
1) Take the average of the voltage multiplied by a reference sine wave at the required frequency, Vs.
2) Take the average of the voltage multiplied by a reference cosine wave at the required frequency, Vc
3) Do the same for current giving Is and Ic.
Now arctan(Vs / Vc) gives the phase of the incoming voltage relative to the internal refence signals.
Similarly arctan(Is / Ic) gives the phase of the incoming current signal.
The difference between the two is the phase shift of current vs voltage.
Your technique will work fine as long as the voltage and current are pure sine waves. There is a problem if there are harmonics, as you mentioned earlier. You can put the current in the imaginary part, but it must be scaled the same as the voltage. Also, I think the resulting phase is half of the actual phase difference, but I might be wrong on that.
Intantaneous power is defined as P = I * V. This needs to be integrated (or summed) over a complete cycle to get the average value for AC power. The power factor is the ratio of the actual power versus the power if I and V are in phase. I think the in-phase power would be the product of the RMS values of V and I.
Dave
Here's the little bit of code I had to write using your awesome FFT object:
repeat timer := cnt repeat i from 0 to 1023 waitcnt(timer += clkfreq/1024) vx[i] := long[voltageSamplePtr]>>1 ix[i] := long[currentSamplePtr]>>1 vy[i] := 0 iy[i] := 0 fft.butterflies(fft#CMD_DECIMATE | fft#CMD_BUTTERFLY, @vx, @vy) fft.butterflies(fft#CMD_DECIMATE | fft#CMD_BUTTERFLY, @ix, @iy) pf := f32.fabs(f32.cos(f32.fsub(f32.atan2(f32.ffloat(vy[60]), f32.ffloat(vx[60])), f32.atan2(f32.ffloat(iy[60]), f32.ffloat(ix[60]))))) long[powerFactorPtr] := pfThanks for the code. One question I had, does it assume a sampling rate of 1024 Hz, or is it just a coincidence that the result locations line up with their frequencies? I've done plenty of FFT in Sysquake and MATLAB, and have never noticed this before.
So what you are saying is that given a bunch of samples of current and voltage, 1024 per cycle each say, we can get the power factor by:
1) Calculating the RMS values of the voltage and current samples, Vrms and Irms.
2) Get RMS power from that: Prms = Vrms * Irms
3) Now calculate the instantaneous power at each sample point Vn * In
4) Take the average of that to get the real power Pr
The Power Factor = Pr / Prms
Further that this is correct even if the waves forms are distorted and full of harmonics.
.
This rather makes the FFT a big hammer to crack a nut:)
Well if you look at it you see that the FFT has no idea about time and no idea about your Hz.
When you see an amplitude in an output bin, bx[60] for example, it just means "60 cycles per 1024 samples". If you happen to be sampling at 1024 samples per second then that represents 60Hz in reality. If you were sampling at 16384 samples per second then bx[60] would represent 960Hz. And so on. I think I have seen this referred to as the "Fourier Frequency" in some literature. It's the frequency within the sample set rather than any real frequency.
Turns out that FFT might be overkill for measuring power factor. But seeing as you have it working nicely I think you should consider calculating the percentage harmonic distortion from the data you have available in the FFT output.
We might guess that the voltage distortion is quite low and constant but the current distortion could be large if there are non-linear loads like rectifiers or SCR controlled light dimmers.
So if you see a high distortion you know it's time to be looking at those harmonics for power factor or consider a different technique.
Actually I would have found that useful to have a few years back when we had trouble getting some equipment working on an UPS. Turned out it was monitoring mains frequency and amplitude for safety reasons but the UPS was generating a lot of harmonics that confused it causing it to trip out all the time.
Anyway, I have now learned more about the mysterious power factor than I ever though I would:)
Your description is much better than the way I was trying to explain it. And, this works with any type of waveform -- not only for pure sinusoids.
Dave
What do you use for a current sense transformer?
RS_Jim
CR8420-1000-G from CR Magnetics. But any should work, as long as you adjust the burden resistance too keep the voltage drop linear with primary current. I'll be posting all the details of the project very soon.
This is more brief than originally intended. Firefox crashed on me 4 times while I was typing.
Not sure if this FFT is what you need but a few comments anyway,
Depends. According to the definition on wikipedia the FFT has no 1/N but the inverse does. Anyway the heater_fft butterfly code does not have 1/N that is left until the magnitude calculation stage.
More importantly the inverse FFT has a different sign for the twiddle factors.
Not a bad idea. I fear the problem here is one of speed. What frequency are you hoping to go up to?
That is a problem. If you want a 1Hz resolution up to 20KHz you are going to need 40K samples per second. What frequency resolution do you think you need?
The "twiddle factor" tables were simply one half cycle of cos and minus sine (wx and wy respectively). They are scaled to +/- 4095. So to double the FFT size we double the size of those tables and recalculate those.
BUT the latest versions of heater_FFT don't define their own sin/cos tables they look up what the need in the Props PROM tables.
My worry is, and I have not worked this through much, that when you double the size of the FFT you also double the max values in the results. Things might start to overflow.
As I said the twiddle factor tables have a range +/-4095. You will see some shifts by 12 in the code where the multiplies are taking care of this fixed point scaling. On could reduce the twiddle
factor range to +/- 2047 or such and reduce the shifts accordingly to reduce overflow problems.
Now if you decide to generate your waves by a simpler route you might like to think about Direct Digital Synthysis(DDS). Here is a chap generating nice sine waves on an AVR with only half a dozen lines of AVR assembler code: http://www.myplace.nu/avr/minidds/index.htm
The audio portion is described on Pages 14 and 15 of his project.
Good thought. The reason I had suggested the methods mentioned above is that using an reverse FT would take enough processing time that the effective sample rate would go through the floor, and while mixing could be accomplished, with only 20MIPs per core I have found that the quality is lost. With the Simpler methods you get just as good of mixing as you would using a reverse FT on a system fast enough to be capable of keeping ahead of the sample rate, with out any loss. and you still have enough time left over to communicate between COGs, thus the sky is the limit. Though I can definitely see showing the wave form on a display, and a reverse FT would be one way to accomplish this as the time constraints are not as important.
My suggestion for using the FFT to display the spectrum of what is being synthesized is just icing on the cake for this project. As you say the time constraints are not there. Just sample, transform and display at whatever speed it goes around at.
P.S. I do like the style of that keyboard in the picture.
With this new code, if you need to calculate magnitude, pass CMD_POWER | CMD_MAGNITUDE. This works exactly same as specifying CMD_MAGNITUDE in original code. If you only pass CMD_POWER, then no square root operation is performed, and on exit bx contains power of each freq component (power is magnitute squared).
Specifying CMD_HAMMING applies the Hamming window to input. If CMD_COMPLEX is also specified, both bx and by are affected by window function, otherwise only bx is affected. Hamming window only works for FFT_SIZE of 8192 and below.
Brilliant, this is a great contribution. I'll roll it into a version 3.0 of heater_fft after I have played with it for a bit.
In case anyone is wondering what this "Hamming Window" is about I think I can summarize it like this:
The FFT works very well if all the frequencies present in the input signal are such that their wavelengths fit exactly into the sample buffer, that is one cycle of a sine wave, say, starting at zero in the first sample, rising to a max, passing through zero in the middle of the buffer, going negative in the second half of the buffer and then returning to zero at the last sample position. Or two, three, four, etc cycles up to the nyquist frequency. With these "neat" waveforms in the input there are no discontinuities in amplitude or slope at the ends of the buffer and a nice clean spectrum is produced by the FFT.
BUT what happens with real world signals where you may have 10.3 cycles in the input? The waves don't fit neatly into the buffer, they don't "wrap around" nicely. There are discontinuities at the ends. These discontinuities appear like noise to the FFT which will now produce an output showing noise frequencies through out the spectrum.
The way to reduce that noise is to reduce the effects of the samples at the ends of the buffers, where the discontinuities occur. This is done by multiplying all the input samples by a function that is zero (or near zero) at the ends of the buffer but rises to a non-zero value at the center of the buffer. There are many such functions. A simple triangular function would work for example. Hamming is such a function based on half a cycle of cosine.
See wikkipedia for more info on this technique.
Now, I have not had much time to check it out yet but I did apply the Hamming window to the test tones that are set up in the heater_fft demo. Thing is with those tones I found when plotting the results that there was not much difference in output. The major difference being that all components are reduced to about 30% amplitude compared to not using the window.
I'm guessing things would look different when looking at lower frequencies and especially fractional frequencies that don't fit as whole cycles into the input buffers.
Andrey, do you have an application in mind for this, I'm sure we'd all like to hear about it.
P.S. This code look really nice in PZST:)
P.P.S. I'm a bit worried about this idea of growing the FFT size passed 1024 samples. There is, I think, a problem that as the FFT size grows the calculations are going to overflow sooner.
I am just playing with various DSP techniques. Using your FFT, I think I could build an audio spectrum analyzer with unique features - 2 independent channels, 31 Hz resolution, virtually any number of frequency channels. And, most exciting - same unit can be used with virtually any display size. Add few shift registers - and make display on 8x8 LED matrix, or 32x16 matrix. Or hook up LCD panel (or even two of them
Just found a stupid bug in the code. Around line 308, it says while it should be
With $2000, a sampled sine wave with applied window function looks like this:
And with $4000:
that is what Hamming window looks like
Before this correction, I was wondering why I see 3rd harmonic bigger than 1st
When you say "coherent signal" do you mean a signal that "fits" nicely in the sample buffer? Sounds like it.
I have also had visions of building such an audio spectrum analyzer or DJ "light organ" from FFT. I just don't have time to keep up with the endless stream of projects that present themselves:) It would be great to see.
There are those who suggest the Hartley Transform is the way to go. Supposedly faster and using less memory. As I don't begin to understand the thing I have not looked at Hartley codes much. It's taken me decades to finally realize how the FFT works:) There is a Hartley Transform for the Prop somewhere.
I'm going to try and create some non-coherent test input for heater_fft that shows off what the Hamming Window can do. Hopefully I can post some plots here soon.
Be aware that the execution time of heater_fft is dependent on input data. That's because the multiply routine has a speed dependent on the size of the numbers being multiplied. Perhaps a fixed time multiply loop would be better.
I'm not sure what I'm looking at there but it does not immediately look right.
I would have thought that a sine wave multiplied by the Hamming function would be alternating equally positive and negative, as it normally does, but with higher amplitude in the middle of the buffer. Your second plot shows a sine wave with a huge DC bias, no negative going.