% Load the frequency domain data
folder='';% = 'C:\Users\haneu\OneDrive\바탕 화면\date\New folder (2)';
filename = '270mvp.csv';
data = readtable(fullfile(folder, filename));
% Extract frequency and FFT amplitude (assumed to be in dB)
f = table2array(data(3:end, 3)); % Frequency data (in Hz)
x_dB = table2array(data(3:end, 4)); % FFT amplitude data in dB
% Number of points in the FFT
N = length(x_dB);
% Compute the sampling frequency (fs)
% Assuming that your frequency axis (f) is spaced evenly
df = f(2) - f(1); % Frequency resolution
fs = N * df; % Sampling frequency
% Convert the amplitude from dB to linear scale
x_linear = 10.^(x_dB / 20);
% Compute the Power Spectral Density (PSD)
pxx = (x_linear.^2) / (fs * N);
% Convert PSD to dB/Hz
pxx_dBHz = 10 * log10(pxx);
% Fundamental frequency
f0 = 1e6; % 1 MHz
% Identify harmonic frequencies
harmonics = f0 * (1:floor(max(f)/f0));
% Remove harmonics from PSD
pxx_dBHz_filtered = pxx_dBHz; % Copy original PSD
for k = 2 %1:length(harmonics)
harmonic_freq = harmonics(k);
[~, idx] = min(abs(f - harmonic_freq)); % Find the closest frequency index
pxx_dBHz_filtered(idx) = -Inf; % Set PSD of harmonics to -Inf (effectively removing)
end
% Plot the original Power Spectral Density (PSD)
figure;
subplot(2, 1, 1);
plot(f/1e6, pxx_dBHz);
xlabel('Frequency (MHz)');
ylabel('PSD (dB/Hz)');
title('Original Power Spectral Density (PSD)');
xlim([0 10])
% Plot the filtered Power Spectral Density (PSD)
subplot(2, 1, 2);
hL=plot(f/1e6, pxx_dBHz_filtered);
xlabel('Frequency (MHz)');
ylabel('PSD (dB/Hz)');
title('Filtered Power Spectral Density (PSD)');
% Highlight fundamental frequency in the plot
hold on;
hL(2)=plot(f0/1e6, pxx_dBHz_filtered(abs(f - f0) < df), 'ro', 'MarkerSize', 8, 'LineWidth', 2);
% Identify harmonic frequencies
harmonics = f0 * (1:floor(max(f)/f0));
% Remove harmonics from PSD
pxx_dBHz_filtered = pxx_dBHz; % Copy original PSD
for k = 2:length(harmonics) % don't do fundamental (first element)
harmonic_freq = harmonics(k);
[~, idx] = min(abs(f - harmonic_freq)); % Find the closest frequency index
pxx_dBHz_filtered(idx) = -Inf; % Set PSD of harmonics to -Inf (effectively removing)
end
% add harmonics visualiztion
xline(harmonics(2:end)/1E6,'r-')
hLg=legend(hL,'Filtered PSD', 'Fundamental Frequency');
xlim([0 10])
You identified the fundamental as well as the harmonics in the harmonics array; when you iterated over the entire array, then you removed it, too. The above to start the loop at two (2) omits the fundamental.
Alternatively, you could define harmonics as
harmonics = f0 * (2:floor(max(f)/f0));
instead and then iterate over the full array from 1:numel(harmonics). That way, the variable name would also reflect the content more exactly than also including the fundamental.
NOTA BENE, however, that the fundamental and the harmonics are not just single-channel peaks; the expanded frequency axis range shows that since the frequency of the fundamental isn't centered exactly at a frequency bin value, the total energy is spread across several bins and so your resulting spectrum still has significant energy content from the harmonics from the spreading of the energy into adjacent bins. The one nearest has been removed, but the neighbors are still there by this "one bin only" technique.
