>>Some comments on your code;=20
>> =20
>> input =3D signal(n) - xn_M;=20
>>
>>This difference between newest and oldest inputs samples is what causes
> the=
>> sinc(x) nulls to be exactly aligned to all the bin frequencies.=20
>>
>>
>> s0 =3D s1 * const_cos - s2 + input;=20
>>
>>This recursion is an iir filter with poles on the unit circle at an
>>angle
> e=
>>qual to the center of the bin frequency.=20
>>
>>Note that at the bin frequency exactly, you get one of the sinc(x) nulls
> al=
>>igned with the pole frequency, so you have zero times infinity, which in
> th=
>>is case gives you a finite number (only practical in the digital
> world!).=
>>=20
>>
>>So what happens if you change the cos term to something in between the
> bins=
>>? Then the problem is that the sinc zeros and the recursive poles don't
> ali=
>>gn anymore, so if you looked at the frequency response of the region
> around=
>> the bin, you would get a peak that was too large (like, infinite), next
> to=
>> a notch that was very close to the peak in frequency.=20
>>
>>One way you MIGHT be able to get what you want is to compute a
>>fractional
> d=
>>elay value around xn_M such that the sinc zeros are shifted so that they
> al=
>>ign with the poles again. But that fractional delay might need to be
>>very
> a=
>>ccurate, and would soak up a fair number if MIPS; there are likely
>>better
> w=
>>ays to tackle the problem. Come to think of it, if you had the ability
>>to
> a=
>>ccurately compute a fractional delay, then you could just use that code
>>as
> =
>>the core of a sample-rate converter, to convert your signal so that the
> des=
>>ired frequency DOES fall into a bin!
>>
>>As others have hinted, if you are trying to estimate the energy of some
> sig=
>>nal in a narrow frequency region, there is nothing "magic" about a
> Goertzal=
>> filter, other than it matches the complex output of a true dft bin. So
>> if
> =
>>you don't really need this mathematically, just use a conventional
> bandpass=
>> filter. Or if you need a complex output, design a complex bandpass
> filter.=
>> In my experience it's rare that you really need to match the dft. It's
> jus=
>>t one of an infinite variety of bandpass filters.=20
>>
>>Bob
>>
> Thank you for your explanation. My idea was that if I know the frequency
> of input signal, I will be able to adapt Goertzel algorithm to compute
> amplitude and phase of single bin (variable, from about 45 to 55 Hz) of
> power line voltage. Samplin frequency is 2400 Hz and length of buffer 48
> samples. So I tried to change angle 2*pi*fp/fs in cos and sin constant
> according frequency fp, but algorithm stopped its function.
> Now I know, why.
> Can you send me tip or link to complex bandpass filter or other
> algorithm, which could be used to compute dft bin instead GA?
>
> Miroslav
Miroslav:
Why didn't you start here? I strongly suggest that you start another
thread, with the topic "Best way to estimate amplitude and phase of
sinusoid".
In the mean time, here's what I suggest, particularly if you can vary
your buffer length slightly according to frequency. This assumes that fp
is known:
Step 1:
Calculate the closest number of bins, from fp: n_bins = 2400Hz / fp.
Step 2:
Calculate the exact frequency according to that number of bins:
fm = 2400Hz / n_bins.
Step 3:
Use some simple algorithm to estimate phase and amplitude from that. I
strongly suggest that you use a sine and cosine wave prototype, i.e., find
x_i = sum(input * cos(th)) / (0.5 * n_bins),
x_q = sum(input * sin(th)) / (0.5 * n_bins)
Now your amplitude is equal to sqrt(x_i^2 + x_q^2), and your phase is
equal to atan2(x_q, x_i).
Your best bet is to compute the angle, th, so that it is zero at the
center of your buffer, and figure that you're computing the phase from
that point in time. This should have the least error due to the
frequency being off.
Your amplitude will be off slightly when your frequency isn't exact, but
it'll be small. It should go as something like sinc(pi * (1 - fn / fm))
or some such -- I don't know what it'll be exactly, so go check. It
should be a small enough error to either ignore, or to calibrate out.
--
Tim Wescott
Wescott Design Services
http://www.wescottdesign.com
Reply by Eric Jacobsen●July 8, 20132013-07-08
On Mon, 08 Jul 2013 16:23:24 -0500, "Mimar" <94571@dsprelated> wrote:
>>Some comments on your code;=20
>> =20
>> input =3D signal(n) - xn_M;=20
>>
>>This difference between newest and oldest inputs samples is what causes
>the=
>> sinc(x) nulls to be exactly aligned to all the bin frequencies.=20
>>
>>
>> s0 =3D s1 * const_cos - s2 + input;=20
>>
>>This recursion is an iir filter with poles on the unit circle at an angle
>e=
>>qual to the center of the bin frequency.=20
>>
>>Note that at the bin frequency exactly, you get one of the sinc(x) nulls
>al=
>>igned with the pole frequency, so you have zero times infinity, which in
>th=
>>is case gives you a finite number (only practical in the digital
>world!).=
>>=20
>>
>>So what happens if you change the cos term to something in between the
>bins=
>>? Then the problem is that the sinc zeros and the recursive poles don't
>ali=
>>gn anymore, so if you looked at the frequency response of the region
>around=
>> the bin, you would get a peak that was too large (like, infinite), next
>to=
>> a notch that was very close to the peak in frequency.=20
>>
>>One way you MIGHT be able to get what you want is to compute a fractional
>d=
>>elay value around xn_M such that the sinc zeros are shifted so that they
>al=
>>ign with the poles again. But that fractional delay might need to be very
>a=
>>ccurate, and would soak up a fair number if MIPS; there are likely better
>w=
>>ays to tackle the problem. Come to think of it, if you had the ability to
>a=
>>ccurately compute a fractional delay, then you could just use that code as
>=
>>the core of a sample-rate converter, to convert your signal so that the
>des=
>>ired frequency DOES fall into a bin!
>>
>>As others have hinted, if you are trying to estimate the energy of some
>sig=
>>nal in a narrow frequency region, there is nothing "magic" about a
>Goertzal=
>> filter, other than it matches the complex output of a true dft bin. So if
>=
>>you don't really need this mathematically, just use a conventional
>bandpass=
>> filter. Or if you need a complex output, design a complex bandpass
>filter.=
>> In my experience it's rare that you really need to match the dft. It's
>jus=
>>t one of an infinite variety of bandpass filters.=20
>>
>>Bob
>>
>Thank you for your explanation. My idea was that if I know the frequency of
>input signal, I will be able to adapt Goertzel algorithm to compute
>amplitude and phase of single bin (variable, from about 45 to 55 Hz) of
>power line voltage. Samplin frequency is 2400 Hz and length of buffer 48
>samples. So I tried to change angle 2*pi*fp/fs in cos and sin constant
>according frequency fp, but algorithm stopped its function.
>Now I know, why.
>Can you send me tip or link to complex bandpass filter or other algorithm,
>which could be used to compute dft bin instead GA?
>
>Miroslav
>Some comments on your code;=20
> =20
> input =3D signal(n) - xn_M;=20
>
>This difference between newest and oldest inputs samples is what causes
the=
> sinc(x) nulls to be exactly aligned to all the bin frequencies.=20
>
>
> s0 =3D s1 * const_cos - s2 + input;=20
>
>This recursion is an iir filter with poles on the unit circle at an angle
e=
>qual to the center of the bin frequency.=20
>
>Note that at the bin frequency exactly, you get one of the sinc(x) nulls
al=
>igned with the pole frequency, so you have zero times infinity, which in
th=
>is case gives you a finite number (only practical in the digital
world!).=
>=20
>
>So what happens if you change the cos term to something in between the
bins=
>? Then the problem is that the sinc zeros and the recursive poles don't
ali=
>gn anymore, so if you looked at the frequency response of the region
around=
> the bin, you would get a peak that was too large (like, infinite), next
to=
> a notch that was very close to the peak in frequency.=20
>
>One way you MIGHT be able to get what you want is to compute a fractional
d=
>elay value around xn_M such that the sinc zeros are shifted so that they
al=
>ign with the poles again. But that fractional delay might need to be very
a=
>ccurate, and would soak up a fair number if MIPS; there are likely better
w=
>ays to tackle the problem. Come to think of it, if you had the ability to
a=
>ccurately compute a fractional delay, then you could just use that code as
=
>the core of a sample-rate converter, to convert your signal so that the
des=
>ired frequency DOES fall into a bin!
>
>As others have hinted, if you are trying to estimate the energy of some
sig=
>nal in a narrow frequency region, there is nothing "magic" about a
Goertzal=
> filter, other than it matches the complex output of a true dft bin. So if
=
>you don't really need this mathematically, just use a conventional
bandpass=
> filter. Or if you need a complex output, design a complex bandpass
filter.=
> In my experience it's rare that you really need to match the dft. It's
jus=
>t one of an infinite variety of bandpass filters.=20
>
>Bob
>
Thank you for your explanation. My idea was that if I know the frequency of
input signal, I will be able to adapt Goertzel algorithm to compute
amplitude and phase of single bin (variable, from about 45 to 55 Hz) of
power line voltage. Samplin frequency is 2400 Hz and length of buffer 48
samples. So I tried to change angle 2*pi*fp/fs in cos and sin constant
according frequency fp, but algorithm stopped its function.
Now I know, why.
Can you send me tip or link to complex bandpass filter or other algorithm,
which could be used to compute dft bin instead GA?
Miroslav
_____________________________
Posted through www.DSPRelated.com
Reply by ●July 4, 20132013-07-04
Some comments on your code;
input = signal(n) - xn_M;
This difference between newest and oldest inputs samples is what causes the sinc(x) nulls to be exactly aligned to all the bin frequencies.
s0 = s1 * const_cos - s2 + input;
This recursion is an iir filter with poles on the unit circle at an angle equal to the center of the bin frequency.
Note that at the bin frequency exactly, you get one of the sinc(x) nulls aligned with the pole frequency, so you have zero times infinity, which in this case gives you a finite number (only practical in the digital world!).
So what happens if you change the cos term to something in between the bins? Then the problem is that the sinc zeros and the recursive poles don't align anymore, so if you looked at the frequency response of the region around the bin, you would get a peak that was too large (like, infinite), next to a notch that was very close to the peak in frequency.
One way you MIGHT be able to get what you want is to compute a fractional delay value around xn_M such that the sinc zeros are shifted so that they align with the poles again. But that fractional delay might need to be very accurate, and would soak up a fair number if MIPS; there are likely better ways to tackle the problem. Come to think of it, if you had the ability to accurately compute a fractional delay, then you could just use that code as the core of a sample-rate converter, to convert your signal so that the desired frequency DOES fall into a bin!
As others have hinted, if you are trying to estimate the energy of some signal in a narrow frequency region, there is nothing "magic" about a Goertzal filter, other than it matches the complex output of a true dft bin. So if you don't really need this mathematically, just use a conventional bandpass filter. Or if you need a complex output, design a complex bandpass filter. In my experience it's rare that you really need to match the dft. It's just one of an infinite variety of bandpass filters.
Bob
Reply by ●July 4, 20132013-07-04
Sme comments on you code;
input = signal(n) - xn_M;
This difference between newest and oldest inputs samples is what causes the sinc(x) nulls to be exactly aligned to all the bin frequencies.
s0 = s1 * const_cos - s2 + input;
This recursion is an iir filter with poles on the unit circle at an angle equal to the center of the bin frequency.
Note that at the bin frequency exactly, you get one of the sinc(x) nulls aligned with the pole frequency, so you have zero times infinity, which in this case gives you a finite number (only possible in the digital world!).
So what happens if you change the cos term to something in between the bins? Then the problem is that the sinc zeros and the recursive poles don't align anymore, so if you looked at the frequency response of the region around the bin, you would get a peak that was too large, next to a notch that was very close to the peak in frequency.
One way you MIGHT be able to get what you want is to do some fractional delay around xn_M such that the sinc zeros are shifted so that they align with the poles again. But that fractional delay might need to be very accurate, and would soak up a fair number if MIPS; there are likely better ways to tackle the problem.
As others have hinted, if you are trying to estimate the energy of some signal in a narrow frequency region, there is nothing "magic" about a Goertzal filter, other than it matches the complex output of a true dft bin. So if you don't really need this mathematically, just use a conventional bandpass filter. Or if you need the complex output, design a complex bandpass filter. In my experience it's rare that you really need to match the dft.
Bob
y_G(n) = (s0 - s1 * const_exp)*exp(-1i*2*pi*k);
s2 = s1; % 1/z
s1 = s0; % 1/z
% y_G(i) = s0 - temp;
xn_M = fifo(1);
end
_____________________________
Posted through www.DSPRelated.com
Reply by Vladimir Vassilevsky●July 4, 20132013-07-04
On 7/3/2013 11:57 AM, Tim Wescott wrote:
> On Wed, 03 Jul 2013 02:09:09 -0500, Mimar wrote:
>
> It's hard to say without knowing the details. Any filter that tries to
> find the amplitude of a sine wave from a non-integer number of cycles is
> going to be subject to difficulties, Goertzel or not.
Huh?
Hint: rotate phasor so it would be integer number of cycles.
Vladimir Vassilevsky
DSP and Mixed Signal Designs
www.abvolt.com
Reply by Tim Wescott●July 4, 20132013-07-04
On Thu, 04 Jul 2013 07:41:11 -0500, Mimar wrote:
>>On Wed, 03 Jul 2013 02:09:09 -0500, Mimar wrote:
>>
>>> Hello somebody,
>>>
>>> does somebody have experiences with sliding Goertzel algorithm? With
>>> K,
>>> which is not integer number? And what about some differences between
>>> G.
>>> filter and G. algorithm?
>>>
>>> Last week I implemented GA in DSP (Freescale 56F8037) and it works
>>> very well for integer K, but for non-integer K this alorithm gives
>>> incorrect results. If I tried to use "normal" GA (without sliding), it
>>> began to work.
>>> Amazing :-). But I cannot understand, why it happens. Any ideas?
>>
>>It's hard to say without knowing the details. Any filter that tries to
>>find the amplitude of a sine wave from a non-integer number of cycles is
>>going to be subject to difficulties, Goertzel or not.
>>
>>Is it possible to post a _succinct_ expression of the algorithm? Not 50
>>lines of code, but maybe 5 lines of math?
>>
>>--
>>
>>Tim Wescott Wescott Design Services http://www.wescottdesign.com
>>
>>Hello,
>
> this a part of my code I have used to simulate of G. algorithm lately.
> It is written in Matlab. So I have already mentioned, for K = integer
> (usually K = 1) it works as well as I want.
>
> %........Goertzel..................................................
> k = 1;
> N = 48: %fs = 2400 Hz pik_term = 2 * pi * k./N;
> const_cos = 2 * cos(pik_term);
> const_exp = cos(pik_term) - 1i*sin(pik_term);
> xn_M = 0;
> fifo = zeros(1, N); % buffer s0 = 0; % states s1 = 0;
> s2 = 0;
> y_G = zeros(1, signal_length);
>
> for n = 1 : 1 : signalľength % loop
>
> fifo = [fifo(2:N), signal(n)]; % sample in buffer
>
> input = signal(n) - xn_M;
> s0 = s1 * const_cos - s2 + input;
> y_G(n) = (s0 - s1 * const_exp)*exp(-1i*2*pi*k);
>
> s2 = s1; % 1/z s1 = s0; % 1/z % y_G(i) = s0 - temp;
> xn_M = fifo(1);
> end
>
> y_G = y_G/N;
>
> y_G = 2*abs(y_G); % result
>
> If I saw this code is functional I overwrote it in C and assembler for
> 56F8300 Freescale core last week. I used fixed point Q7.24 format to get
> good accuracy. But the problem with real number K continues still :-(.
Well, that's hardly a succinct algorithm description.
It appears that you're implementing the transfer function
H(z) = z^2 / (z^2 - 2 * cos * th * z + 1)
over however many samples takes your fancy, and saving the result into a
vector. You don't say what you're doing with the vector, so we can't
comment on that part of your algorithm.
So, by not explicitly using any window at all you're using a rectangular
window. What are you doing with y_G? The most simple-minded things to
do with it would make your algorithm sensitive to both the phase of the
input as well as the period not being an integer.
You really need to cough up more information. And, since the Goertzel
algorithm _isn't_ the most efficient way to do a one-bin DFT on a chip
that has a good MAC instruction, perhaps you could share with us what you
really want to get out of the exercise, and then we could help you reach
your real goal efficiently.
--
Tim Wescott
Wescott Design Services
http://www.wescottdesign.com
Reply by Mimar●July 4, 20132013-07-04
>On Wed, 03 Jul 2013 02:09:09 -0500, Mimar wrote:
>
>> Hello somebody,
>>
>> does somebody have experiences with sliding Goertzel algorithm? With K,
>> which is not integer number? And what about some differences between G.
>> filter and G. algorithm?
>>
>> Last week I implemented GA in DSP (Freescale 56F8037) and it works very
>> well for integer K, but for non-integer K this alorithm gives incorrect
>> results. If I tried to use "normal" GA (without sliding), it began to
>> work.
>> Amazing :-). But I cannot understand, why it happens. Any ideas?
>
>It's hard to say without knowing the details. Any filter that tries to
>find the amplitude of a sine wave from a non-integer number of cycles is
>going to be subject to difficulties, Goertzel or not.
>
>Is it possible to post a _succinct_ expression of the algorithm? Not 50
>lines of code, but maybe 5 lines of math?
>
>--
>
>Tim Wescott
>Wescott Design Services
>http://www.wescottdesign.com
>
>Hello,
this a part of my code I have used to simulate of G. algorithm lately. It
is written in Matlab. So I have already mentioned, for K = integer (usually
K = 1) it works as well as I want.
%........Goertzel..................................................
k = 1;
N = 48: %fs = 2400 Hz
pik_term = 2 * pi * k./N;
const_cos = 2 * cos(pik_term);
const_exp = cos(pik_term) - 1i*sin(pik_term);
xn_M = 0;
fifo = zeros(1, N); % buffer
s0 = 0; % states
s1 = 0;
s2 = 0;
y_G = zeros(1, signal_length);
for n = 1 : 1 : signalľength % loop
fifo = [fifo(2:N), signal(n)]; % sample in buffer
input = signal(n) - xn_M;
s0 = s1 * const_cos - s2 + input;
y_G(n) = (s0 - s1 * const_exp)*exp(-1i*2*pi*k);
s2 = s1; % 1/z
s1 = s0; % 1/z
% y_G(i) = s0 - temp;
xn_M = fifo(1);
end
y_G = y_G/N;
y_G = 2*abs(y_G); % result
If I saw this code is functional I overwrote it in C and assembler for
56F8300 Freescale core last week. I used fixed point Q7.24 format to get
good accuracy. But the problem with real number K continues still :-(.
_____________________________
Posted through www.DSPRelated.com
Reply by dszabo●July 3, 20132013-07-03
>Hello somebody,
>
>does somebody have experiences with sliding Goertzel algorithm? With K,
>which is not integer number? And what about some differences between G.
>filter and G. algorithm?
>
>Last week I implemented GA in DSP (Freescale 56F8037) and it works very
>well for integer K, but for non-integer K this alorithm gives incorrect
>results. If I tried to use "normal" GA (without sliding), it began to
work.
>Amazing :-). But I cannot understand, why it happens. Any ideas?
>
>Thank you a lot.
>
>Miroslav
>
Could you elaborate on what you mean by incorrect results?
_____________________________
Posted through www.DSPRelated.com
Reply by robert bristow-johnson●July 3, 20132013-07-03
On 7/3/13 11:26 AM, Tim Wescott wrote:
> On Wed, 03 Jul 2013 11:12:58 -0700, robert bristow-johnson wrote:
>
>> On 7/3/13 9:57 AM, Tim Wescott wrote:
>>> On Wed, 03 Jul 2013 02:09:09 -0500, Mimar wrote:
>>>
>>>
>>> It's hard to say without knowing the details. Any filter that tries to
>>> find the amplitude of a sine wave from a non-integer number of cycles
>>> is going to be subject to difficulties, Goertzel or not.
>>>
>>>
>> isn't that what window functions are for?
>
> Well, one question to ask is if he's using a window function.
he's gotta be using *some* window function. if he doesn't know, then
it's likely the rect(). and a sliding rect() window applied to a
sinusoid that does not have an integer number of cycles in the width of
that rect() "is going to be subject to difficulties". it will matter
much in how that rect() window lands on top of the waveform.
but with a decent window with decaying tails on both sides, it might not
matter so much how the window lands on top of the waveform. and if he's
tracking a known frequency (that is not at the center of the DFT bin),
he should be able to use the sliding DFT for the closest DFT bin and the
adjacent bins and interpolate results. but then, if he knew the exact
frequency, the dot product that he's doing with the DFT could be using
the exact frequency, rather than the nice DFT bin frequencies.
> That's why I was asking for details!
details are for the anal retentive.
:-)
--
r b-j rbj@audioimagination.com
"Imagination is more important than knowledge."