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.

  • $\begingroup$ The article seems to be behind a paywall... $\endgroup$ – Ben Feb 16 at 12:45
  • $\begingroup$ 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. $\endgroup$ – Ben Feb 16 at 12:51
  • $\begingroup$ Thanks, then How should I generate the H and do convolution? $\endgroup$ – ASiDesigner Feb 16 at 14:33
  • $\begingroup$ you can check out the Paper in here sci-hub.do/10.1109/tpwrs.2019.2892599 $\endgroup$ – ASiDesigner Feb 16 at 15:00
  • $\begingroup$ The filter is a moving-average, h[n] is explained on page 2 $h_{av}$ $\endgroup$ – 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;.

  • $\begingroup$ 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? $\endgroup$ – ASiDesigner Feb 18 at 11:16
  • $\begingroup$ No, the rest of the processing required is explained in the paper. You need to shift the frequency before filtering. $\endgroup$ – Ben Feb 18 at 12:34
  • $\begingroup$ 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 $\endgroup$ – ASiDesigner Feb 20 at 6:54
  • $\begingroup$ How should I do it in matlab? $\endgroup$ – ASiDesigner Feb 23 at 7:55
  • $\begingroup$ Which part specifically? $\endgroup$ – 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

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

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


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.