# A Fast Power Grid Frequency Estimation Approach Using Frequency Shift Filtering

I have read the paper A Fast Power Grid Frequency Estimation Approach Using Frequency Shift Filtering.

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,

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


Now How should I calculate the hp(n) values using MATLAB? Finally I should convolve the result with this sample sin wave to get 50.1 Hz estimation.

• The article seems to be behind a paywall... – Ben Feb 16 at 12:45
• From the look of it, they seem to perform the Park transform with a fixed-frequency (50 Hz) in your case. Since, the frequency is not exactly 50 Hz, the d-q signals won't be DC but will have a frequency of f - 50 Hz or 0.1 Hz in your case. – Ben Feb 16 at 12:51
• Thanks, then How should I generate the H and do convolution? – ASiDesigner Feb 16 at 14:33
• you can check out the Paper in here sci-hub.do/10.1109/tpwrs.2019.2892599 – ASiDesigner Feb 16 at 15:00
• The filter is a moving-average, h[n] is explained on page 2 $h_{av}$ – Ben Feb 16 at 16:54

Basically they cascade moving-average filters, or CAF.

a Moving-average filter has the following coefficients

$$h_{av}(m) = 1/M, m = 0, 1.. M-1$$

In Matlab it would look like this

h = ones(1,M)./M;
y = filter(h,1,x);


If you want to cascade three moving-average filters

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


M should be equal to the number of samples per period. So if your frequency is approximately 50 Hz and your sampling frequency is 12800 Hz then M = 256;.

Edit : I modified your script to make it work. There were 2 problems, first of all you only processed 256 samples, your couldn't measure anything significant because the transient response of your 3 cascaded filters is 768 samples long. Secondly, you filtered the reference and not your complex signal...

clear
clc
%ADC sin wave, it represent the data that I have captured with the ADC
%and want to estimate the frequency
grid on
hold on

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

%calculating  xs

%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)

• Thanks, for the info, so how hp is calculated, so we do the convolution and estimate that this example X signal has 50.1Hz frequency? – ASiDesigner Feb 18 at 11:16
• No, the rest of the processing required is explained in the paper. You need to shift the frequency before filtering. – Ben Feb 18 at 12:34
• Dear Ben thanks for the feedback, would you please tell me How to do it in matlab? I'm a bit confused how to do it. thanks – ASiDesigner Feb 20 at 6:54
• How should I do it in matlab? – ASiDesigner Feb 23 at 7:55
• Which part specifically? – Ben Feb 23 at 14:21

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

clear
clc
%ADC sin wave, it represent the data that I have captured with the ADC
%and want to estimate the frequency
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

%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?

• You're doing many things wrong. For starters you filter only 256 samples with FIR samples that have 256 coefficients. You will always be in the transient response. Secondly, why do you filter the reference 3 times? You need to filter xs with 3 cascaded moving-average filters. I modified your script and it works, I get a frequency of 50.1 Hz. – Ben Apr 29 at 11:46
• Dear ben would you please share your code. – ASiDesigner May 2 at 8:02