Why am i getting huge values on low frequencies of my PSD ?

Hi,

I'm using an IMU on a stepper motor which provide rotations a 1Hz. I retrieve pose quaternion as well as timestamp via processing. Then, I calculate a rotation matrix from the quaternion and retrieve the angle of my sensor on the horizontal plane according to a reference vector. This signal is theorically periodic with a frequency of 1Hz. Then, I perform FFT and PSD analysis to get the fundamental frequency of my signal on the angle of rotation. I also try Welch analysis.

However, my PSD value contains weird values on low frequencies. I don't know why i'm getting these values. I suppose it's due to the stepper motor which generate micro vibrations due to gears but i want to be sure that my algorithm is correct.

Here is my code :

   if true
 %import csv file
 M = csvread('3.csv');

 qx = M(:,1);
 qy = M(:,2);
 qz = M(:,3);
 qw = M(:,4);
 t = M(:,5);

 %size of M (same size for each column)
 s = size(qx, 1);
 %because of loop, need s/2
 s2 = fix(s/2);

 %get the angle from quaternion
 %create orthogonal vector of Z
 axis = [1, 0, 0];
 ang = zeros(s,1);
 for i = 1:s
     qxx = qx(i);
     qyy = qy(i);
     qzz = qz(i);
     qww = qw(i);
     %create rotation matrix from quaternion
     a00 = 1 - 2 * qyy * qyy - 2 * qzz * qzz;
     a01 = 2 * (qxx * qyy - qww * qzz);
     a02 = 2 * (qxx * qzz + qww * qyy);
     a10 = 2 * (qxx * qyy + qww * qzz);
     a11 = 1 - 2 * qxx * qxx - 2 * qzz * qzz;
     a12 = 2 * (qyy * qzz - qww * qxx);
     a20 = 2 * (qxx * qzz - qww * qyy);
     a21 = 2 * (qyy * qzz + qww * qxx);
     a22 = 1 - 2 * qxx * qxx - 2 * qyy * qyy;
     %rotation matrix from world frame to module frame
     A = [a00, a01, a02; a10, a11, a12; a20, a21, a22];
     %apply rotation matrix to this vector
     %(coordinates expressed in the world frame)
     % A\b = inv(A)*b
     rotated = A\axis';
     %project rotated vector to the horizontal plane
     flattened = [rotated(1), rotated(2), 0];
     flattened = flattened/norm(flattened);
     %get the angle
     ang(i,1) = acos(dot(axis, flattened));
 end 

 figure
 plot(ang);
 % computes the fast fourier transform of M for angle
 Ygz = fft(ang);

 %Compute the power spectral density for angle
 PSDYgz = Ygz.*conj(Ygz)/s;

 %acquisition frequency
 hz = 30;
 %calculate x axis (frequency)
 f = hz/s*(1:s2);

 %find peaks of the frequency sample for angle
  [pksgz,locsgz] = findpeaks(PSDYgz(1:s2), 'SortStr', 'descend', 'NPeaks', 1);

 figure  

 hold on
 plot(f, PSDYgz(1:s2))
 plot(f(locsgz), pksgz, 'Marker' , 'v')
 str = sprintf('%0.3e', f(locsgz));
 str2 = strcat(num2str(str), ' Hz');
 h = text(f(locsgz), -0.1 ,str2);
 set(h, 'HorizontalAlignment', 'center', 'VerticalAlignment', 'top');
 str = sprintf('%0.3e', pksgz);
 str2 = strcat(num2str(str), ' g²/Hz');
 h = text(-0.5, pksgz, str2);
 set(h, 'HorizontalAlignment', 'right'); 

 title('Power Spcetral Density for gyroscope')
 xlabel('Frequency (Hz)')
 ylabel('Power (g²/Hz)')
 hold off

 %perform welch algorithm
 [pxx, ff] = pwelch(ang, hz);
 %plot welch
 figure
 plot(ff, 10*log10(pxx));

   end

Angle plot : good at the beginning and weird at the end

PSD plot : big values on low frequencies

Welch plot : big values on low frequencies

Thanks a lot,

Best.

Answer 1

[Mona Mahboob Kanafi]

Hello,

We don't know what you expect to see in your power spectrum. You just mentioned you see some weird high peaks in your PSD. This is quite normal as the peak in zero frequency is related to the mean of your signal which is not zero. Second, your welch plot in in logarithmic scale, in order to see your true PSD check it also in logarithmic scale. Third, check the number of points in low frequency, probably you don't have that many peaks in low frequency, the first one is at zero frequency related to the mean and then two more maybe related to general low frequency pattern of the input signal. Fouth, if your signal is not periodic in nature, you must apply some window function to your signal, before applying FFT, due to spectral leakage:

e.g. http://www.ni.com/white-paper/4844/en/

[Maxence]

So, the peak in 0 frequency is due to the mean of my signal. If I divide all my value by the mean of my signal, I should get rid of this first peak right ?

Then, why need I plot the PSD on logarithmic scale ?

And you're right, I do not have many values on low frequencies. I'm trying to apply a FIR high pass filter at 0.1Hz on my signal to get rid of the first peaks. I suppose I need to apply the filter on my input signal values and not on the FFT values right ? I know these peaks are due to low frequencies but i would like to know the true reason of these peaks. Is a sensor error or motor low freqency vibration? May be it's due to the gyroscope drift or the gyroscope inertia ?

Finally, i'm gonna read your link to try the window function. Thank you very much.

[Mona Mahboob Kanafi]

First, to remove the mean:

ang = ang - mean(ang); % now the mean is zero

You need to plot your PSD in logarithmic scale, because from your pwelch results which are actually in logarithmic scale, you can see the peaks in high frequencies. To make it more clear, if you take a look at your signal, what FFT does is to show you the multiple waves (and their corresponding frequencies) which are added together to construct your signal. Based on your signal, you can see that the amplitude of your high frequency components are quite small compared to low frequency behaviour; thus, you could detect them in logarithmic scale, as the components have orders of magnitude difference in energy.

You need to apply filter in frequency domain, not time domain. Filter design is an other story, but as an example refer to this link:

http://www.mathworks.com/matlabcentral/answers/9647-use-fft-x-as-a-highpass-filter

About the nature of the peaks you get from the experiments, I can't be of help as this is not area of my expertise.

Hope these help,

-Mona

[Maxence]

I applied the mean on my signal and of course it removed the 0 Hz peak. That's a good point.

However, even if I apply my high pass filter (0.1Hz cutoff) on frequency domain, I have always my low frequencies peaks which are not correct. I have a peak jsut on 1Hz which the good one. Plot on logarithmic scale reduce the peak amplitude difference but does not solve the problem :)

I'm not sure but, can my initial condition on the angle integration be a problem ? I set my first angle to 0 but I think it's not 0. May be FFT can consider this value as a very low frequency wave ?

Thanks for all.

[Soares Guiamar]

I had the same problem. Try to use a rectangular window in pwelch(pwelch(x,rectwin(length(x)),....). For me, it worked.

[Oskar Kilgus]

Hey @Maxence, im facing quite a similar problem and i´m wondering if you found a proper solution back in 2020.

I tried getting rid of the offset and also tried a rectwin as @Soares Guiamar suggested, but still there are these low-frequency peaks/the roll off from 0 Hz.

Glad for any reply, thanks!