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

75 views (last 30 days)
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.

Answers (1)

Mona Mahboob Kanafi
Mona Mahboob Kanafi on 12 Feb 2016
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:
  5 Comments
Soares Guiamar
Soares Guiamar on 17 Jul 2020
I had the same problem. Try to use a rectangular window in pwelch(pwelch(x,rectwin(length(x)),....). For me, it worked.
Oskar Kilgus
Oskar Kilgus on 24 May 2023
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!

Sign in to comment.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!