how to use matrix pencil or prony in matlab

Hello,

I have a distorted frequency data with respect to time in an excel file. I can graph it out in Matlab.

Since it is distorted frequency, I want to find how many external frequencies are combined with the fundamental frequency of 50Hz. After that I want to find the damping value with respect to frequency.

I want to use Matrix Pencil or Prony's method to solve this problem but I don't know where to start. Can someone guide me?

Answer 1

[Milan Bansal]

Hi uzzi

Here is how you can use the Prony's method to design a filter to find the external frequency and damping values.

1. Import your data from the Excel file into MATLAB.
2. Ensure your data is centered around zero if necessary and remove any trends or noise.
3. Apply Matrix Pencil or Prony's Method. Both methods are used for estimating the parameters of a sum of damped sinusoids.
4. Check how many frequencies are close to the fundamental frequency of 50 Hz and identify additional frequencies.
5. Analyze the damping values corresponding to each frequency.
6. Plot the original and reconstructed signals to verify the decomposition.

Here is a basic implementation using Prony's method.

% Load data from Excel
[data, txt, raw] = xlsread('filename.xlsx');

% Assume the first column is time and the second column is frequency
time = data(:, 1);
frequency = data(:, 2);

% Detrend the data
frequency = detrend(frequency);

% Define the order of the model (number of damped sinusoids)
modelOrder = 10; % Adjust based on your data

% Use Prony's method
[a, g] = prony(frequency, modelOrder, modelOrder);

% Calculate poles
poles = roots(a);

% Calculate frequencies and damping factors
frequencies = abs(angle(poles)) / (2 * pi * mean(diff(time)));
damping = -log(abs(poles)) / mean(diff(time));

% Display results
disp('Frequencies:'); disp(frequencies);
disp('Damping values:'); disp(damping);

% Reconstruct the signal
reconstructed_signal = zeros(size(frequency));
for k = 1:length(poles)
    reconstructed_signal = reconstructed_signal + g(k) * exp(poles(k) * time);
end

% Plot
figure;
subplot(2, 1, 1);
plot(time, frequency);
title('Original Signal');
ylabel('Frequency');

subplot(2, 1, 2);
plot(time, real(reconstructed_signal));
title('Reconstructed Signal');
ylabel('Frequency');

Please refer to the following documentation links to learn more:

• detrend: https://www.mathworks.com/help/matlab/ref/detrend.html
• prony: https://www.mathworks.com/help/signal/ref/prony.html

Hope this helps!