Frequency spectral analysis, how do I calculate it?

Hi guys,

I am currently trying to run a frequency analysis on acceleration signal, to look at the harmonic ratio.

To work out the harmonic ratio, the sum of even harmonics are divided by sum of odd harmonics.

The code below is I modified from an example of MATLAB. So what I am trying to do is, trying to add the value of peaks at each frequency(e.g. 1Hz, 2Hz). However, as my data is around 2257 sample points with 200 FS. Thus, my frequency bin is 200/2257.

It would be helpful if I can get some advice on how to get harmonics from the signal.

THanks!

% Sampling frequency = 200;

x_temp = input_data;

L = length(x_temp); % Length of signal

Nfft = 2^(nextpow2(L)); % The length of the signal should be a power of 2

f2 = abs(fft(x_temp,Nfft)); % two-sioded spectrum f2

f1(2:L/2) = 2*f2(2:L/2);

power_f1 = f1.^2;

f = (0:L/2-1)*(fs/L);

plot(f,power_f1)

xlabel('Frequency')

ylabel('Power')

% Edit: loading figure as this is easier than downloading to view.
openfig('FFt_1.fig');

Answer 1

[William Rose]

@Yong Kim,

The signal in the file you uploaded seems to be truncated. See the plot below.

The idea of harmonics only makes sense when there is a clearly known or obs=rvable fundamental frequency in th signal. You r signal does not have an ovbvios fundamental frequency. Maybe the fundamental frequency would be more obvious if the signal were not truncated. See the result of the code below. I renamed your .mat file so it would not conflict with other file names. I included the renamed file in my answer so that the code would be runnable in this window.

%freqAnalysisMLA20220206.m

clear
load('freqAnalysisMLA20220206.mat');  %file provided by user
% Sampling frequency 
fs=200;
x = input_data;
L = length(x);          %Length of signal  
t=(0:L-1)/fs;           %time vector
window=L; noverlap=0; nfft=L;
[Pxx,f]=cpsd(x,x,window,noverlap,nfft);
figure;
subplot(211), plot(t,x,'-rx');
xlabel('Time (s)'); ylabel('x(t)'); grid on
subplot(212), plot(f,Pxx,'-b.')
xlabel('Frequency (Hz)'), ylabel('Power'), grid on; xlim([0,0.5])

I have used cpsd() from the signal processing toolbox. Let me know if you do not have access to this toolbox. The spikes in the power spectrum suggest a fundamental frequency of f0=0.1058/4 (frequency of 4th peak=0.1058, divided by four).

Notice that the peaks do not exactly line up with the frequency spacing. The first peak has its power evenly split between two frequencies, while the power at some ofther peaks is concentrated in a single frequency sample. Thisleads you to ask: How shall I compute the spectrum at each harmonic? Anwer 1: Take the value of the spectrum at the frequency closest to the harmonic. Answer 2: Integrate the area under the curve for a certain range of frequencies near the harmonic. If you go with Answer 2, you must decide how many points to use and again, what to do about the frequency spacing.

I have not limited the signal to a power of two number of points, because there is no good reason to do it. It used to be done for speed reasons, but computers are so fast now, and these data sets are so tiny, that speed is not an issue. When you zero-pad or truncate a signal to reach a power of two, you alter the sampling in the frequency domain, for no good reason. If you pad with zeros, the spectrum has a higher apparent resolution than is really justified by the data. If you truncate, you are throwing away useful information.

[William Rose]

@Yong Kim,

Te estimate derived above for funamental frequency was wrong, because the truncaiton of the data created spikes in the spectrum corresponding to the recpirocal of the duraiton of the untruncated part of the signal.

The signal was all zeros from point 956 onward, as shown in the plot of input_data, below. I have highlighted the points which are peaks, because they show successive repeats of the signal.

The time (in samples) of the peaks are

tpeak=[159,402,635,870];

The times between peaks are

disp(diff(tpeak))

   243   233   235

The average duration of one cycle is

avgPer=(tpeak(4)-tpeak(1))/3;

The average period, in seconds, is

fs=200; Tavg=avgPer/fs; disp(Tavg)

    1.1850

and therefor an estimate of the fundamental frequency is

f0=1/Tavg; disp(f0);

    0.8439

Please upload an un-truncated version of the signal.

[William Rose]

I mde a mistake. I forgot to pass fs whn calling cpsd(). I did

[Pxx,f]=cpsd(x,x,window,noverlap,nfft);

and I should have done

[Pxx,f]=cpsd(x,x,window,noverlap,nfft,fs);

This caused the values on my frequency axis to be incorrect. After correcting for this, the frequencies in the power spectrum make more sense, as shown below.

clear
load('freqAnalysisMLA20220206.mat');  %file provided by user
% Sampling frequency 
fs=200;
x = input_data;
L = length(x);          %Length of signal  
t=(0:L-1)/fs;           %time vector
window=L; noverlap=0; nfft=L;
[Pxx,f]=cpsd(x,x,window,noverlap,nfft,fs);
figure;
subplot(211), plot(t,x,'-rx');
xlabel('Time (s)'); ylabel('x(t)'); grid on
subplot(212), plot(f,Pxx,'-b.')
xlabel('Frequency (Hz)'), ylabel('Power'), grid on; xlim([0 10])

My estimate of the fundamental frequency, from the plot above, is

f0=3.347/4; disp(f0)  % freq. of 4th peak, divided by 4

    0.8367

This is very close to the estimate f0=0.8439 Hz, which I gave above, and which was based on the time-domain appearance of the signal.

Sorry for the error.

Answer 2

[William Rose]

@Yong Kim,

This code computes "harmonic ratio, the sum of even harmonics are divided by sum of odd harmonics". I define even harmonics as multiples 2,4,6,8,10 of the fundamental, and odd harmonics as multiples 1,3,5,7,9 of the fundamental. I stopped after 10 because there is very little pwer beyond the 10th harmonic. The code also plots the harmonics 1 through 10. See console output and plot below.

%freqAnalysisMLA20220206.m
%Yong Kim with edits from W. Rose, 2022-02-06
clear
load('freqAnalysisMLA20220206.mat');  %file provided by user
fs=200;                             % sampling frequency 
%x = input_data;
x = input_data(1:956);              % non-zero part of input_data
L = length(x);                      % length of signal  
t=(0:L-1)/fs;                       % time vector
window=L; noverlap=0; nfft=L;
[Pxx,f]=cpsd(x,x,window,noverlap,nfft,fs);
figure;
subplot(211), plot(t,x,'-rx');      %time-domain signal
xlabel('Time (s)'); ylabel('x(t)'); grid on
subplot(212), plot(f,Pxx,'-b')     %power spectrum
xlabel('Frequency (Hz)'), ylabel('Power'), grid on; xlim([0,20]);
hold on;
[pks,locs]=findpeaks(Pxx,'MinPeakHeight',0.1);
f0=f(locs(4))/4; 
fprintf('Fundamental frequency=%.3f Hz.\n',f0);

Fundamental frequency=0.837 Hz.

feven=(2:2:10)*f0;              %first 5 even harmonic frequencies
fodd=(1:2:9)*f0;                %fundamental and 4 odd harmonic frequencies
oddfreqindex=zeros(1,5);
evenfreqindex=zeros(1,5);
for i=1:5
    [~,oddfreqindex(i)]=min(abs(fodd(i)-f));
    [~,evenfreqindex(i)]=min(abs(feven(i)-f));
end
plot(f(oddfreqindex),Pxx(oddfreqindex),'gs','LineWidth',2);
plot(f(evenfreqindex),Pxx(evenfreqindex),'r^','LineWidth',2);
legend('Power','Odd Harmonics','Even Harmonics');

PowerOdd=sum(Pxx(oddfreqindex));
PowerEven=sum(Pxx(evenfreqindex));
HarmonicRatio=PowerEven/PowerOdd;
fprintf('Harmonic Ratio=Even/Odd Harmonic Power=%.3f.\n',HarmonicRatio)

Harmonic Ratio=Even/Odd Harmonic Power=0.588.

Try it. Good luck with your work.

[William Rose]

@Yong Kim,

Some of my code may have been hard to understand. Therefore I am providing comments below for some of the lines. The method I used to compute f0, the fundamental frequency, could have found a frequency that was not exactly one of the frequencies in the power spectrum. Therefore I wrote code that would find the index of each frequency in the spectrum that is closest to each computed harmonic frequency. This is why the code I posted above is a bit complicated. I hope the comments I added below will explain this.

%next: find peaks in spectrum
[pks,locs]=findpeaks(Pxx,'MinPeakHeight',0.1);
%next: fundamental freq=frequency of 4th peak divided by 4
f0=f(locs(4))/4;
fprintf('Fundamental frequency=%.3f Hz.\n',f0);
feven=(2:2:10)*f0;         %first 5 even harmonic frequencies
fodd=(1:2:9)*f0;           %fundamental and 4 odd harmonic frequencies
oddfreqindex=zeros(1,5);   %initialize array to contain indices of harmonics
evenfreqindex=zeros(1,5);  %initialize array to contain indices of harmonics
for i=1:5
    %next: find the indices of the frequencies in the spectrum that are closest
    %to the calculated odd and even harmonic frequencies
    [~,oddfreqindex(i)]=min(abs(fodd(i)-f));
    [~,evenfreqindex(i)]=min(abs(feven(i)-f));
end
%next: plot symbols at the even and odd harmonics of the spectrum
plot(f(oddfreqindex),Pxx(oddfreqindex),'gs','LineWidth',2);
plot(f(evenfreqindex),Pxx(evenfreqindex),'r^','LineWidth',2);
legend('Power','Odd Harmonics','Even Harmonics');
PowerOdd=sum(Pxx(oddfreqindex));    %total power at harmonics 1,3,5,7,9
PowerEven=sum(Pxx(evenfreqindex));  %total power at harmonics 2,4,6,8,10
HarmonicRatio=PowerEven/PowerOdd;   %ratio: even/odd power