- Start with N pairs (x1, y1), (x2, y2), …, (xN, yN).
- Compute the means (xm, ym).
- Compute the covariance matrix Cov= (cxx , cxy ; cyx , cyy) where
- cxx= [ ∑(xi-xm)² ] / (N-1)
- cxy = cyx = [ ∑(xi-xm)(yi-ym) ] / (N-1)
- cyy = [ ∑(yi-ym)² ] / (N-1)
- Find the eigenvalues (λ1, λ2) of the covariance matrix (λ1 > λ2).
- Find the unit eigenvectors (v11 ; v21) and (v12 ; v22) of the covariance matrix.
- The best fit ellipse has its center at (xm, ym).
- The semimajor axis of the 95% ellipse has direction (v11; v21) and length sqrt(5.99 λ1).
- The semiminor axis of the 95% ellipse has direction (v12; v22) and length sqrt(5.99 λ2).
- The area of the 95% ellipse is A = 5.99 π sqrt(λ1 λ2)
- The angle between the semimajor axis and the y-axis is θ = sin-1(v11).

# Principal axis + phase and inclination.

10 views (last 30 days)

Show older comments

Hello, I would appreciate any help/assistance in calculating the following variables using the attached code.:

Semi-major axis (Mayor), Semi-minor axis (Minor), Inclination and Phase

of Horizontal velocities u and v are represented as complex numbers complex(u,v) with the purpose of plotting the ellipse as shown in the attached figure.

%horizontal velocity ellipse rotating clockwise.

% Example data for time series of u and v velocities

u = [ -3.4500 -2.3500 0.8500 -1.3500 0.2000 0.7500 3.8500 5.0000 4.9000 8.0000 6.2000 14.6500 21.8000 14.8500 13.0500];

v = [ 13.1000 12.3000 8.0000 3.8000 2.9000 1.8500 3.7500 1.0000 3.3000 4.8500 -3.9000 9.9000 11.4000 16.5500 20.3000];

% Combine u and v velocities into complex vector

wind_complex = complex(u, v);

% Compute covariance matrix

cov_matrix = cov(real(wind_complex), imag(wind_complex));

% Find eigenvalues and eigenvectors

[eigenvectors, eigenvalues] = eig(cov_matrix);

% Major and minor axes

major_axis_length = max(diag(eigenvalues));

minor_axis_length = min(diag(eigenvalues));

% Calculate inclination and phase

inclination = atan2(eigenvectors(2,1), eigenvectors(1,1)); % angle in radians

phase = atan2(mean(imag(wind_complex)), mean(real(wind_complex))); % angle in radians

% Convert inclination and phase to degrees

inclination_degrees = rad2deg(inclination);

phase_degrees = rad2deg(phase);

% Display results

disp(['Major Axis Length: ', num2str(major_axis_length)]);

disp(['Minor Axis Length: ', num2str(minor_axis_length)]);

disp(['Inclination (degrees): ', num2str(inclination_degrees)]);

disp(['Phase (degrees): ', num2str(phase_degrees)]);

Major Axis Length: 67.5505

Minor Axis Length: 29.7218

Inclination (degrees): -54.0886

Phase (degrees): 51.446

>>

##### 0 Comments

### Answers (1)

William Rose
on 22 Apr 2024

I assume you want to do a best-fit ellipse to the wind data. Here is a link to my notes on how to find the best-fit ellipse to a set of x,y points. The notes refer to center-of-presssure data , but it could be any set of x,y points.

I apologize for the poor quality of the file. This old file online is the only copy of my notes that I can locate at the moment. Key points from the notes:

The following procedure finds the “best fit” ellipse.

Plot your wind data:

% Example data for time series of u and v velocities

u = [-3.45 -2.35 0.85 -1.35 0.20 0.75 3.85 5.00 4.90 8.00 6.20 14.65 21.80 14.85 13.05];

v = [13.10 12.30 8.00 3.80 2.90 1.85 3.75 1.00 3.30 4.85 -3.90 9.90 11.40 16.55 20.30];

figure; plot(u,v,'k*');

grid on; xlabel('u (m/s)'); ylabel('v (m/s)'); hold on

Calculate best fit ellipse. I notice that you do not take the square root of the eigenvalues, in your calculations. You should. And you should scale them by a suitable factor, such as sqrt(5.99) to catch 95% of the points, or sqrt(4.61) to catch 90% of the points. (The scaling factors come from the chi squared distribution, and are based on the assumption that the noise in the data is Gaussian and independent for u and v.)

##### 1 Comment

William Rose
on 22 Apr 2024

Here is a script that computes the semimajor and semiminor axis length and the inclination of the majopr axis with respect to the +x axis. It makes the figure below.

The script also displays these results:

>> bestFitEllipse

Major Axis Length = 20.14.

Minor Axis Length = 13.36.

Inclination = -144.1 degrees.

I am not sure what you mean by the phase angle. If the data lay on or close to an ellipse, and if there were an associated time vector, then I suppose the phase could be defined as the angle of the vector from the ellipse center to the point at time zero, or to the first point in the data set. But this data do not appear to lie on or close to an ellipse.

### See Also

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!