如何利用MATLAB中Fourier变换去除小时风速数据的周期性

Got it, let's break down exactly how to use Fourier series to strip out those daily, weekly, monthly, and yearly periodic patterns from your 7-year hourly wind speed data in MATLAB, then build your prediction model. I'll walk you through each step with actionable code examples:

First, get your data into MATLAB and convert timestamps into a continuous hour index—this makes modeling periodicity over time much simpler. Assuming your data is stored in a CSV with datetime and wind speed columns:

% Load your dataset (replace with your actual file path)
wind_data = readtable('hourly_wind_speed.csv');

% Convert datetime to a continuous hour count since the first measurement
start_time = wind_data.DateTime(1);
hours_since_start = hours(wind_data.DateTime - start_time);

% Extract the wind speed values
wind_speed = wind_data.WindSpeed;

Map the daily, weekly, monthly, and yearly cycles to hours (matching your hourly data frequency):

• Daily: 24 hours
• Weekly: 24*7 = 168 hours
• Monthly: ~730.49 hours (average, accounting for varying month lengths)
• Yearly: ~8766 hours (average, including leap years)

Fourier series models cyclic patterns as sums of sine and cosine terms. We'll create a design matrix of these terms, then use least squares to fit the periodic signal. You can adjust the number of harmonics to capture more complex cyclic behavior:

% Define periods in hours
periods = [24, 168, 730.49, 8766];
n_harmonics = 2; % Start with 2 harmonics per period; adjust later if needed

% Build the design matrix for Fourier terms
design_matrix = ones(length(hours_since_start), 1); % Constant offset term
for T = periods
    for k = 1:n_harmonics
        % Add cosine and sine terms for each harmonic
        design_matrix = [design_matrix, ...
            cos(2*pi*k*hours_since_start/T), ...
            sin(2*pi*k*hours_since_start/T)];
    end
end

% Fit the periodic model using least squares
coeffs = design_matrix \ wind_speed;
periodic_signal = design_matrix * coeffs;

% Remove the periodic component from the original data
detrended_wind = wind_speed - periodic_signal;

Plot the original data, extracted periodic signal, and detrended data to confirm cyclic patterns are properly removed:

figure('Position', [100, 100, 800, 600]);
subplot(3,1,1);
plot(wind_data.DateTime, wind_speed);
title('Original Hourly Wind Speed');
ylabel('Wind Speed');
grid on;

subplot(3,1,2);
plot(wind_data.DateTime, periodic_signal);
title('Extracted Periodic Components (Day/Week/Month/Year)');
ylabel('Periodic Signal');
grid on;

subplot(3,1,3);
plot(wind_data.DateTime, detrended_wind);
title('Wind Speed After Removing Periodicity');
ylabel('Detrended Wind Speed');
grid on;

With detrended data in hand, you can use any time series prediction model. Here's an example using ARIMA (a popular choice for wind speed forecasting):

% Convert detrended data to a timeseries object
ts = timeseries(detrended_wind, wind_data.DateTime);

% Fit an ARIMA model (adjust p,d,q based on your data's autocorrelation)
arima_model = arima(2, 1, 1); % AR(2), Integrated(1), MA(1)
fit_model = estimate(arima_model, ts.Data);

% Forecast the next 24 hours of wind speed
[forecast_vals, ~, ~] = forecast(fit_model, 24);

% Plot the forecast against recent historical data
figure;
plot(ts.Time(end-100:end), ts.Data(end-100:end), 'b-', ...
    ts.Time(end) + hours(1:24), forecast_vals, 'r--');
legend('Recent Detrended Data', '24-Hour Forecast');
title('Wind Speed Prediction (Post-Periodicity Removal)');
ylabel('Wind Speed');
grid on;

• Tune Harmonics: If the periodic signal doesn't capture all cyclic behavior, increase n_harmonics (try 3 or 4).
• Exact Periods: For more accuracy, calculate exact monthly/yearly periods using MATLAB's datetime functions instead of averages (handles leap years and varying month lengths).
• Check Stationarity: Use the adftest function to confirm your detrended data is stationary—this is critical for models like ARIMA.
• Alternative Decomposition: If Fourier series isn't working well, try MATLAB's built-in decompose function for seasonal-trend decomposition.

内容的提问来源于stack exchange，提问作者L N