Fast Power Grid Frequency Estimation

[ali_asadzadeh]

Hi,
I have read the attached paper.
I want to prove it using a numerical example, I have problem how to calculate the hp(n) values,
Let's suppose that I have a sample signal which I have sampled with an ADC and for simplicity I will use Matlab to make the signal. (the signal has a 50.1HZ frequency) and I want to estimate it with this paper approach,

Code: [Select]

fs = 12800;                    % Sampling frequency (samples per second)
dt = 1/fs;                   % seconds per sample
StopTime = 0.02;             % 20 mili seconds or a 1 full cycle of 50Hz AC
t = (0:dt:StopTime-dt)';     % seconds
F = 50.1;                    % Sine wave frequency (hertz)
X = sin(2*pi*F*t);           % Generate sine wave

acording to paper it seems that I should cascade some moveing avrage filters like this code

Code: [Select]

M = 256;
h = ones(1,M)./M;
y1 = filter(h,1,X);
y2 = filter(h,1,y1);
y3 = filter(h,1,y2);

Then I calculated the convolution of the X and y3 to determine the frequency

Code: [Select]

F = conv2(X, y3);

And now getting frequency from two samples of the F, like this

Code: [Select]

DF = ((F(1) - F(256)) * 256 * 50 )/ (2*pi*255)

So estimated f should be like this

Code: [Select]

Fest = 50 + DF;
Which is 24.0013 and it’s far from the 50.1 answer according to the paper results.

So the question is that am I doing it in the right way? Because I do not get some parts of the paper.

Any help would be highly appreciated

[ali_asadzadeh]

I need some help, please help me find what am I doing wrong?

[ali_asadzadeh]

I may find some hints about why I'm getting wrong results, the question is how to generate a complex sin wave in matlab?

[rstofer]

Here's a simple example of a fundamental along with the 3rd and 5th harmonic.

Code: [Select]

  t = 0:0.001:2*pi;
  y = 0.636620 * sin(t) + 0.212207 * sin(3*t) + 0.17324 * sin(5*t);
  plot(t,y)
[/font]

It approximates a square wave.  The coefficients have been scaled for some reason I fail to recall.  Nevertheless, the shape is right.

[ali_asadzadeh]

Thanks rstofer, do you have any Idea what's wrong with my code in matlab? How should I evaluate the results of the paper in matlab?

[ali_asadzadeh]

So, Now I have done a little bit more, Now I'm generating a complex reference sine wave, The things that I have done are as follow

Code: [Select]

clear
clc
%ADC sin wave, it represent the data that I have captured with the ADC
%and want to estimate the frequency
ADCsin = dsp.SineWave(1,50.1);
ADCsin.ComplexOutput = false;
ADCsin.SampleRate  = 12800;
ADCsin.SamplesPerFrame = 256;
ADC = ADCsin();
plot(ADC);
grid on
hold on

%ref sin wave with 50Hz complex value
Refsin = dsp.SineWave(1,50);
Refsin.ComplexOutput = true;
Refsin.SampleRate  = 12800;
Refsin.SamplesPerFrame = 256;
ref = Refsin();
plot(imag(ref))

%calculating  xs
xs = ref .* ADC;

%calculate the CAF
M = 256;
h = ones(1,M)./M;
y1 = filter(h,1,ref);
y2 = filter(h,1,y1);
y3 = filter(h,1,y2);

%doing the actual final convolution
F = conv2(xs, y3);

%calculating the final result
Fest = 50 - (angle(F(4)) - angle(F(3)) * 256*50)/(2*pi*90)

the answer is clearly wrong! so what I have done wrong?

[penfold]

Is your application of the filter correct? Do you get the same result as formally defining it by repeated convolutions of h_avg as the paper does, theres a part of me saying that the resulting length of the filter affects the result.. do you get a similar result from a single application of the filter() function as you do with two applications or is it entirely different?

[ali_asadzadeh]

I get different answers, you can test the code in Matlab and give feedback if you find something, thanks.

[nctnico]

I'm not quite sure whether this is such a good method. From a quick glance: What it basically does is demodulating the frequency offset into a near DC offset which is prone to being influenced by harmonics without being able to filter these properly.

Regarding the wrong result: I'm missing the arg() from the final calculation in your code. I have not dug too deep into the paper but is seems it is using the angles to determine the frequency offset.

[ali_asadzadeh]

I'm not quite sure whether this is such a good method. From a quick glance: What it basically does is demodulating the frequency offset into a near DC offset which is prone to being influenced by harmonics without being able to filter these properly.

Regarding the wrong result: I'm missing the arg() from the final calculation in your code. I have not dug too deep into the paper but is seems it is using the angles to determine the frequency offset.

Isn't the arg means the angle? or am I missing something?

[penfold]

I get different answers, you can test the code in Matlab and give feedback if you find something, thanks.

I'm not running Matlab (I'm only on gnu octave) so your revised code doesn't play well without me installing dsp toolbox etc

When you say different answers, do you mean different answers from convolved vs successive filter applications or different for each filter order.

Have you tried using the graph of F = conv2(xs, y3) to see how its affecting the results? I noticed there sign change with the original (incorrect) code which suggests that there's something very wrong with the filter application.

[ali_asadzadeh]

I'm not running Matlab (I'm only on gnu octave) so your revised code doesn't play well without me installing dsp toolbox etc

When you say different answers, do you mean different answers from convolved vs successive filter applications or different for each filter order.

Have you tried using the graph of F = conv2(xs, y3) to see how its affecting the results? I noticed there sign change with the original (incorrect) code which suggests that there's something very wrong with the filter application.

If you have an octave code that works, Maybe I can test it hear in MATLAB too. and I mean from different results if I use y1,y2, or y3 for convolution, besides the paper says it has about mHz accuracy and none of the answers are even close

[nctnico]

I'm not quite sure whether this is such a good method. From a quick glance: What it basically does is demodulating the frequency offset into a near DC offset which is prone to being influenced by harmonics without being able to filter these properly.

Regarding the wrong result: I'm missing the arg() from the final calculation in your code. I have not dug too deep into the paper but is seems it is using the angles to determine the frequency offset.

Isn't the arg means the angle? or am I missing something?

You are right. I overlooked that. Angle() should do the same as arg() (Why didn't they call it arg() ?)

I took a look at the paper again. I think you need to carefully go through each step and check whether it is right. Note that the paper isn't the final version so it may contain errors or is missing a step. I think the CAF filter may need more research on how this should be implemented. An easier path could be to find the final version or a paper which describes the same method in a more clear way.

[gf]

Try this:

Code: [Select]

  M = 256;
  N = 256*M;
  fs = 50*M;
  t = [0:N-1]/fs;
  d_t = t(2);
  % reference signal
  ref = exp(-1i*2*pi*50*t);
  % frequency sweep 50-52 Hz
  f = 50 + 2*[0:N-1]/N;
  % input signal
  sig = sin(2*pi*f.*t);
  % down-convert
  bbsig = sig.*ref;
  % 3rd oder moving average filter
  bbsig = conv(bbsig, ones(1,M)/M);
  bbsig = conv(bbsig, ones(1,M)/M);
  bbsig = conv(bbsig, ones(1,M)/M);
  % remove settling of filter
  bbsig = bbsig(4*M:end-4*M);
  d_angle = diff(angle(bbsig));
  % unwrap angle wrap-around from 2*pi -> 0
  d_angle(d_angle<-pi) += 2*pi;
  d_angle(d_angle>pi) -= 2*pi;
  % still unclear why this is off by a factor of 2 ???
  % d_f = d_angle/(2*pi*d_t);
  d_f = d_angle/(2*pi*d_t) / 2;
  plot(50+d_f);
  ylim([45 55])
  title("frequency vs time")
  grid on

IMO the key is that the algorithm is supposed to operate on a cotinuous stream of samples, not on a single block of size M. I simulate this stream by a much larger block and discard leading and trailing samples after filtering which are garbage until the filter has settled. The calculated frequency difference is still off by a factor of 2. I did not find the mistake yet . So I added an extra division by 2. Apparentyl that works, but that's just a hack (I'm not happy if things work with a hack, and I don't understand why).

[ali_asadzadeh]

Thanks for the code gf,
this part of code does not work in matlab

Code: [Select]

d_angle(d_angle<-pi) += 2*pi;
  d_angle(d_angle>pi) -= 2*pi;what did you want to do? so I can redoe it in matlab

[gf]

Thanks for the code gf,
this part of code does not work in matlab

Code: [Select]

d_angle(d_angle<-pi) += 2*pi;
  d_angle(d_angle>pi) -= 2*pi;what did you want to do? so I can redoe it in matlab

The code works in Octave. I guess the operators += and -= are possibly not understood by Matlab. They have the same meaning as in C, i.e. A+=B is a short form for A=A+B.

[gf]

The calculated frequency difference is still off by a factor of 2. I did not find the mistake yet . So I added an extra division by 2. Apparentyl that works, but that's just a hack (I'm not happy if things work with a hack, and I don't understand why).

Did some more investigations. The previously assumed need for an extra division by 2 turned turned out to be definitively wrong (which is a good sign - it would not fit with my phi=onega*t worldview). So I removed it again. For constant frequencies the result looks actually correct now. It rather seems to be the transient response which looks pretty strange (see plot, orange is actual frequency, blue is measured frequency). After the transition, it returns to the correct value, though. Strange

[gf]

Found the problem now. My simulated input signal

Code: [Select]

sig = sin(2*pi*f.*t);
contains phase discontinuities. The method does not like them, but they are unlikely in a real power grid. Replacing it with a correct, phase-continuous chirp

Code: [Select]

f(1) = 0; sig = sin(cumsum(2*pi*f*d_t));
solves the problem and the output looks like expected now
Actually I find this an interesting way to implement a kind of streaming hilbert transformer

EDIT: Attached is a plot of the (simulated) measured frequency, starting at 51Hz, stepping down do 49Hz, ramping up to 51Hz again, and staying at 51Hz.
Risetime of step response is about 500 samples, i.e. approx. two 50Hz-cycles or ~40ms.

Here is the latest version of my Octave code:

Code: [Select]

  M = 256;
  N = 256*M;
  fs = 50*M;
  t = [0:N-1]/fs;
  d_t = t(2);
  % reference signal
  ref = exp(-1i*2*pi*50*t);
  % simulated frequency profile:
  % start at 51Hz, step down to 49Hz, ramp up to 51Hz
  % and remain at 51Hz
  f = [repmat(51,1,N/4) repmat(49,1,N/4) repmat(51,1,N/2)];
  f(N/2:N/2+16384-1) = 49 + 2*[0:16384-1]/16384;
  % create phase-continuous chirp signal for the frequency profile
  phi = 2*pi*f.*d_t; phi(1) = 0;
  sig = sin(cumsum(phi));
  % down-convert
  bbsig = sig.*ref;
  % 3rd oder moving average filter
  bbsig = conv(bbsig, ones(1,M)/M);
  bbsig = conv(bbsig, ones(1,M)/M);
  bbsig = conv(bbsig, ones(1,M)/M);
  % discard settling of filter
  bbsig = bbsig(3*M:end-3*M);
  d_angle = diff(angle(bbsig));
  % unwrap angle wrap-around from 2*pi -> 0
  d_angle(d_angle<-pi) += 2*pi;
  d_angle(d_angle>pi) -= 2*pi;
  % calculate instantaneous frequency
  d_f = d_angle/(2*pi*d_t);
  plot(50+d_f);
  ylim([45 55])
  title("frequency vs time")
  grid on


[ali_asadzadeh]

Dear gf
Hi,
Thanks for the feedback ,would you please resend your new complete code, also would you please use a Complete code for += and -= so it can be done in matlab too, so I can test it here too.

[gf]

The complete code of my simulation was already attached to the previous message. It is not more than that. A streaming DSP implementation would of course look different, processing one sample at a time when it arrives, and not operating on blocked data. The Octave assignment operators are documented here. I leave the trivial conversion to Matlab code as an exercise for the "Digital Expert from 8-bits to 64-bits" . Octave syntax is very similar to Matlab anyway. In order to understand the algorithm described in the paper you need to understand the concepts of Analtic Signals and Instantaneous Phase and Frequency, how to generate an analytic signal from a real signal, and a few DSP basics. One property of the chosen SMA filter is that is has zeros at N*50 Hz, cancelling harmonics at N*50Hz (and still suppressing N*(50Hz+df) pretty well, for small df). With a filter order of only 2, there is still noticable ripple from the -(100Hz+df) mixing product left in the filtered signal, so order 2 is not enough. With an order of 3-4, suppression is already pretty good, up to frequency deviations of a few Hz. You could even go higher, but drawback of a high order is that it reduces the response time - so it's a trade-off.

[ali_asadzadeh]

thanks, I'm a little bit confused, can we do it on a 256 samples only, where is the estimated frequency in your code? my final goal is to fully understand this approach and try it on an STM32

[gf]

The algorithm per se does not require any particular blocks, but can be reazied as a DSP which works on an continuous stream of samples, captured at a fixed sampling rate. For each sample you retrieve from the ADC and feed into the DSP engine, you can get one frequency estimate out.

OTOH, the length of each of the 3-4 cascaded moving average filter stages must indeed be 20ms (i.e. exactly one period of the nominal mains frequency). For a sampling rate of 12800 Sa/s, each filter stage needs to have 256 taps then. The condition also implies that the sampling rate must be an integral multiple of the nominal mains frequency, too.

The simulation is not implemented as a streaming DSP, but rather in an "Octave/Matlab-friendly" way. It creates 64k samples (~5 seconds of) of simulated input data with a given, simulated frequency profile (in the variable f), and applies the steps of the algorithm to the whole vector in parallel (as it is usual practice in Matlab/Octave). The simulated input signal is the vector sig, and the vector of estimated frequencies is 50+d_f, which is eventually plotted at the end (-> one estimated frequency per sample).

[ali_asadzadeh]

I managed to edit the matlab code, now it's working as expected,

Code: [Select]

clear
clc
%ADC sin wave, it represent the data that I have captured with the ADC
%and want to estimate the frequency
ADCsin = dsp.SineWave(1,50.124);
ADCsin.ComplexOutput = false;
ADCsin.SampleRate  = 12800;
ADCsin.SamplesPerFrame = 768;
ADC = ADCsin();
plot(ADC);
grid on
hold on

%ref sin wave with 50Hz complex value
Refsin = dsp.SineWave(1,50);
Refsin.ComplexOutput = true;
Refsin.SampleRate  = 12800;
Refsin.SamplesPerFrame = 768;
ref = Refsin();
plot(imag(ref))

%calculating  xs
xs = ref .* ADC;

%calculate the CAF
M = 256;
h = ones(1,M)./M;
y1 = filter(h,1,xs);
y2 = filter(h,1,y1);
y3 = filter(h,1,y2);


%calculating the final result
Fest = 50 - (angle(y3(end)) - angle(y3(end-1))) * 12800/(2*pi)

Now I want to know if it's possible to make the algorithm simpler, is there a way to combine all the 3 filters to a pre-calculated fixed stream to only do the convolution one time,
I mean is there a way to combine all these 3 lines of code to a one time filter?

Code: [Select]

%calculate the CAF
M = 256;
h = ones(1,M)./M;
y1 = filter(h,1,xs);
y2 = filter(h,1,y1);
y3 = filter(h,1,y2);
So that we do the convolution only one time?

[gf]

Now I want to know if it's possible to make the algorithm simpler, is there a way to combine all the 3 filters to a pre-calculated fixed stream to only do the convolution one time,
I mean is there a way to combine all these 3 lines of code to a one time filter?

Code: [Select]

%calculate the CAF
M = 256;
h = ones(1,M)./M;
y1 = filter(h,1,xs);
y2 = filter(h,1,y1);
y3 = filter(h,1,y2);
So that we do the convolution only one time?

The combinedd filter kernel is h3=conv(h,conv(h,h)).
But there is no computational saving at the end, since h3 is (almost) 3x as long as h then.

For a µC implementation, a runnung average filter has the advantage that you can calculate it with only one addition and one subtraction per sample (and storage space for keeping the last M samples). For three filters in series, it is still only 3 additions and 3 subtractions per sample then. OTOH, convolution with h3 were (3*M-2) multiplications and (3*M-2) additions per sample, i.e. much more expensive, particularly if M is large.

[ali_asadzadeh]

Thanks gf, is there some clever way of combining the first convolution into the 3 stage filter, and get a vector so, all the work could be done with a single run of conv?

I mean combining the ref signal with the 3 stage filters.