FFT is shifting my artificial signal?

I'm working on a function that performs an FFT on a signal that comes from a weather station (The primary function of which is to separate a pressure signal into its mean, wave, and perturbation elements. This is refered to as a Reynold's Triple Decomposition).

The function is complete, but I need to validate the function. So I created an artificial signal that is comprised of a mean component, 3 dominant oscillations, and finally a random variable that will represent our perturbations. The artificial signal is generated using this code:

%% Inputs 
n = 18000;          % Number of values in signal
m = 2;              % Range of random element 
f1 = 0.01;          % Frequency of primary resonance (Hz)
f2 = 0.2; 
f3 = 3; 
i1 = 0.5;           % Intensity of wave components
i2 = 0.5; 
i3 = 0.5; 
avgVal = 500;       % Mean value for signal
%% Setting up signal
X = 1:1:n; 
Y_raw = ones(1,n); 
Y_mean = Y_raw + (avgVal-1); 
Y_wave = ((Y_raw .* (i1.*sin(f1.*X))) + (Y_raw .* (i2.*sin(f2.*X))) + (Y_raw .* (i3.*sin(f3.*X))));
Y_MpW = Y_wave + Y_mean; 
Y_rand = -m + (m + m)*rand([1,n]);
Y_final = Y_MpW + Y_rand;

Then I have a second function that loads the data from the first as if it were a normal input.

%% Importing/Organizing data 
load('Y_final.mat'); 
m_1 = length(Y_final);
%% Extraction of mean parameter (Step 1)
AvgElement = mean(Y_final); 
DataMinusMean = Y_final - AvgElement;
%% Creation of the Hanning window 
H_per = 0.05; 
H_bound = m_1 * H_per; 
H_NoTop = hann(H_bound*2); 
H_Top = ones(1, m_1); 
H_Top(1, 1:H_bound) = H_NoTop(1:H_bound); 
H_Top(1, (m_1-H_bound):m_1) = H_NoTop(H_bound:end); 
%% Period of the primary wave (Step 2)
Fs = 20;              % Sampling frequency    
nf = Fs/2;            % Nyquist Frequency
T = 1/Fs;             % Sampling period       
L = m_1;              % Length of signal
t = (0:L-1)*T;        % Time vector
f = Fs*(0:(L/2))/L;
FT_P_dcs = fft(H_Top.*DataMinusMean);
Fu_P_dcs = ((FT_P_dcs./length(DataMinusMean)).*conj(FT_P_dcs./length(DataMinusMean)));
% conversion to the energy spectrum
n = (Fs.*(1:nf))./L;
DeltaN = Fs/L; 
Ea_dcs = 2*Fu_P_dcs; 
Spp_dcs = Ea_dcs./DeltaN;  
P2_dcs = abs(Spp_dcs/L); 
P1_dcs = P2_dcs(1:L/2+1); 
P1_dcs(2:end-1) = 2*P1_dcs(2:end-1); 
%% PLotting stuff
figure(1); 
semilogx(f,(P1_dcs)) 
hold on; 
xline(0.01, 'r-', 'LineWidth', 2);
xline(0.2, 'r-', 'LineWidth', 2);
xline(3, 'r-', 'LineWidth', 2); 
grid on; 
title("DCS", 'FontSize', 20);
xlabel("f (Hz)", 'FontSize', 20)
ylabel("Intensity", 'FontSize', 20)
hold off; 

Now, the problem....

I have artificially generated a signal that has 3 distinct sinusoidal functions. One at 0.01 Hz, one at 0.2 Hz, and one at 3 Hz. When I plot my FFT (Along with 3 vertical lines at those corresponding points) it is very clear that there is something wrong. The FFT is picking up on 3 very strong frequencies but they are all about 3 times the expected frequency.

0.01 -> 0.3222

0.2 -> 0.6367

3.0 -> 9.5489

Here is an image of the plot:

Any sugestions for how I could go about fixing this?

Answer 1

[Paul]

Hi Zach,

a) create the X vector with spacing based on the sampling frequency, Fs

b) create the wave vector with sinusoids computed from 2*pi*f1, etc. because f1 is in Hz, not rad/sec

%% Inputs 
n = 18000;          % Number of values in signal
m = 2;              % Range of random element 
f1 = 0.01;          % Frequency of primary resonance (Hz)
f2 = 0.2; 
f3 = 3; 
i1 = 0.5;           % Intensity of wave components
i2 = 0.5; 
i3 = 0.5; 
avgVal = 500;       % Mean value for signal
%% Setting up signal

%X = 1:1:n;
Fs = 20;
X = (0:n-1)/Fs;

Y_raw = ones(1,n); 
Y_mean = Y_raw + (avgVal-1);

%Y_wave = ((Y_raw .* (i1.*sin(f1.*X))) + (Y_raw .* (i2.*sin(f2.*X))) + (Y_raw .* (i3.*sin(f3.*X))));
Y_wave = ((Y_raw .* (i1.*sin(2*pi*f1.*X))) + (Y_raw .* (i2.*sin(2*pi*f2.*X))) + (Y_raw .* (i3.*sin(2*pi*f3.*X))));


Y_MpW = Y_wave + Y_mean; 
Y_rand = -m + (m + m)*rand([1,n]);
Y_final = Y_MpW + Y_rand;


m_1 = length(Y_final);
%% Extraction of mean parameter (Step 1)
AvgElement = mean(Y_final); 
DataMinusMean = Y_final - AvgElement;
%% Creation of the Hanning window 
H_per = 0.05; 
H_bound = m_1 * H_per; 
H_NoTop = hann(H_bound*2); 
H_Top = ones(1, m_1); 
H_Top(1, 1:H_bound) = H_NoTop(1:H_bound); 
H_Top(1, (m_1-H_bound):m_1) = H_NoTop(H_bound:end); 
%% Period of the primary wave (Step 2)
Fs = 20;              % Sampling frequency    
nf = Fs/2;            % Nyquist Frequency
T = 1/Fs;             % Sampling period       
L = m_1;              % Length of signal
t = (0:L-1)*T;        % Time vector
f = Fs*(0:(L/2))/L;
FT_P_dcs = fft(H_Top.*DataMinusMean);
Fu_P_dcs = ((FT_P_dcs./length(DataMinusMean)).*conj(FT_P_dcs./length(DataMinusMean)));
% conversion to the energy spectrum
n = (Fs.*(1:nf))./L;
DeltaN = Fs/L; 
Ea_dcs = 2*Fu_P_dcs; 
Spp_dcs = Ea_dcs./DeltaN;  
P2_dcs = abs(Spp_dcs/L); 
P1_dcs = P2_dcs(1:L/2+1); 
P1_dcs(2:end-1) = 2*P1_dcs(2:end-1); 
%% PLotting stuff
figure(1); 
semilogx(f,(P1_dcs)) 
hold on; 
%xline(0.01, 'r-', 'LineWidth', 2);
%xline(0.2, 'r-', 'LineWidth', 2);
%xline(3, 'r-', 'LineWidth', 2); 
grid on; 
title("DCS", 'FontSize', 20);
xlabel("f (Hz)", 'FontSize', 20)
ylabel("Intensity", 'FontSize', 20)
hold off; 

[Zach Ruble]

Of course it's a simple math error.

Thanks Paul!