FFT is shifting my artificial signal?
4 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
Zach Ruble
am 27 Sep. 2023
Kommentiert: Zach Ruble
am 27 Sep. 2023
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:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1494252/image.png)
Any sugestions for how I could go about fixing this?
0 Kommentare
Akzeptierte Antwort
Paul
am 27 Sep. 2023
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;
Weitere Antworten (0)
Siehe auch
Kategorien
Mehr zu Discrete Fourier and Cosine Transforms finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!