[SOLVED] measuring sine wave distotion in matlab

[Lekshmi B S]

Is there any method for calculating the distortion of a sine wave,, in matlab....?
When i did extended kalman filter estimation with input as sine wave, i got the predicted and updated output sine waves as distorted..and i would like to find the amount of distortion to limit the output...So please help if anyone has any idea about this.....

Thank you in advance....

 

[ArNvC]

Use FFT in MATLAB, if signal frequency is 'f', then 2*f, 3*f ...up to N*f inside the frequency band (or inband) are harmonic. Now, sum the total power of harmonic, then you get total harmonic power / signal power is harmonic distortion in db.

Or, use formula : % THD = sqrt(harmonic power)/sqrt(signal power)

 

[Lekshmi B S]

Hi arunava7,

This is the code i am using..While plotting apriori_state and aposteriori_state, i am getting sine waves which are distorted..The distortion of that waveform needs to be calculated..I dont have any other information about the waves,, naa..Then how can i find the amount of distortion/change in shape to the particular signals.....

N = 350;
W = [ 1 0 ; 0 1 ]; % State-error
V = [ 1 0 ; 1 0 ]; % Measurement-error
H = [ 1 0 ; 0 0 ]; % Measurement-state
Q = [ 0.001 0 ; 0 0 ]; % Process noise covariance
R = [ 0.1 0 ; 0 0.01 ]; % Measurement noise covariance

x = zeros(2,N);
apriori_state = zeros(2,N);
aposteriori_state = zeros(2,N);
apriori_error = zeros(2,2,N);
aposteriori_error = zeros(2,2,N);
z = zeros(2,N);
K = zeros(2,2,N);

x(:,1) = [ 1 ; 3*pi/500 ];
aposteriori_state(:,1) = [ 1 ; 1*pi/500 ];
aposteriori_error(:,:,1) = [ 1 0 ; 0 1 ];
randn('state',2);
n1=sqrt(Q)*[randn;randn];
n2=sqrt(R)*[randn;randn];

for i = 2:N    
    x(:,i) = [sin(x(2,i-1)*(i-1));x(2,i-1)] + n1;
    z(:,i) = x(:,i) + n2;
    apriori_state(:,i) = [sin(aposteriori_state(2,i-1)*(i-1));aposteriori_state(2,i-1)];
    Fi = [0 (i-1)*cos(aposteriori_state(2,i-1)*(i-1)) ; 0 1 ];
    Qi = Q;
    Ri = R;
    apriori_error(:,:,i) = Fi*aposteriori_error(:,:,i-1)*Fi' + Qi;
    K(:,:,i) = apriori_error(:,:,i)*H' / (H*apriori_error(:,:,i)*H'+Ri);
    residual(:,i) = (z(:,i) - apriori_state(:,i));
    aposteriori_state(:,i) = apriori_state(:,i) + K(:,:,i) * residual(:,i);
    aposteriori_error(:,:,i) = (eye(2) - K(:,:,i)*H) * apriori_error(:,:,i);
    error1(:,i)=apriori_state(:,i)-x(:,i);
    error2(:,i)=aposteriori_state(:,i)-x(:,i);
end

 

[Lekshmi B S]

Also how should i find the frequency and power of these sine waves to calculate the THD in the way as you suggested..?
Please do reply.....It is somewhat urgent.....

Thank you....

 

[ArNvC]

You do have the equation of the sine wave? then you apply fourier transform.

 

[Lekshmi B S]

Hi arunava7,
Can you please check whether this code will be enough to find the THD using your suggestion..?

XX1 = abs(fft(aposteriori_state)).^2;

% to find the maximum signal power
XX2 = sort(XX1);
power = find(XX1(1:700)==XX2(700));
signal_power = XX1(power)  


harmonic_power = sum([2:N])/350  % to find the sum of harmonic power

THD_percentage= sqrt(harmonic_power)/sqrt(signal_power)

By using this code, i got the value for THD as 43.9.. It is in percentage,, naa..? So is it enough?... And what is the normal expected range for the THD of the sine wave?...........
Thank you...

 

[ArNvC]

Ya, you have used the formula for %THD. It is difficult to say the normal range for THD for a sine wave but it is a clean sinewave with THD usually below 3% in households.

Anyway, if you can measure the total power (usually this is done in time domain) and the dc and 1st (funtamental) components, you have
Harmonics_power = Total_power - DC_power - Fundamental_power. This is another method but it is error prone if harmonics are less.

 

[Lekshmi B S]

Ok thank you arunava7 for your help...