find 3rd harmonics of the signal in matlab

% Generate a signal with 3rd harmonics
Fs = 1000;              % Sampling frequency
t = 0:1/Fs:1;           % Time vector
f1 = 50;                % Fundamental frequency
f3 = 3*f1;              % 3rd harmonic frequency
s = sin(2*pi*f1*t) + sin(2*pi*f3*t);    % Signal with 3rd harmonics

% Compute the Fourier transform of the signal
n = length(s);          % Length of signal
Y = fft(s);             % Compute Fourier transform
Pyy = Y.*conj(Y)/n;     % Power spectral density

% Find the frequency index of the 3rd harmonic
f = Fs*(0:(n/2))/n;     % Frequency vector
f3_idx = find(f == f3);

% Plot the power spectral density and mark the 3rd harmonic
figure
plot(f,Pyy(1:n/2+1))
hold on
plot(f(f3_idx), Pyy(f3_idx), 'ro')
title('Power spectral density')
xlabel('Frequency (Hz)')
ylabel('Power/Frequency (dB/Hz)')
816 chars
25 lines