Continuing from the last article, we have obtained a collection result: an array named allData. This array has a size of [4096×180], containing all the wind speed data collected over 180 seconds. The entire process can be found on the following webpage
- Using Matlab to Operate Hanghua CTA04 Hot-Wire Anemometer 1: Parameter Settings
- Using Matlab to Operate Hanghua CTA04 Hot-Wire Anemometer 2: Automatic Calibration
- Using Matlab to Operate Hanghua CTA04 Hot-Wire Anemometer 3: Manual Calibration
- Using Matlab to Operate Hanghua CTA04 Hot-Wire Anemometer 4: Data Collection
You can refer to our series of public articles from 2023 to process the collected data
- Hot-Wire Data Analysis Method 1: Mean Speed and Turbulence Intensity
- Hot-Wire Data Analysis Method 2: Frequency Spectrum Analysis
- Hot-Wire Data Analysis Method 3: Wavenumber Spectrum Analysis
- Hot-Wire Data Analysis Method 4: Autocorrelation and Integral Scale
- Hot-Wire Data Analysis Method 5: Minimum Turbulence Scale
This article summarizes the programs from the above five articles for your reference. The content includes calculations of mean, frequency spectrum, wavenumber spectrum, autocorrelation function, integral scale, Taylor scale, Kolmogorov scale, etc. This series is complete. The next series will involve using Python to operate the Hanghua CTA04 anemometer. Stay tuned.
%Author:N. Gao, [email protected], Hanghua Inc www.hanghualab.comclear all;close all;
Sample_Frequency=4096;Sample_Number_Per_Block=4096;Velocity_Loaded=allData;t=[1:Sample_Number_Per_Block]/Sample_Frequency; % Time corresponding to each sample point in the data block
T=Sample_Number_Per_Block/Sample_Frequency;df=1/T;f=[0:Sample_Number_Per_Block-1]*df;Number_Blocks=180;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Block_Selected=40; % Randomly select the 40th group from these 180 groups Velocity_Instantaneous=Velocity_Loaded(:,Block_Selected); figure, plot(t,Velocity_Instantaneous); xlabel('t, s');ylabel('u, m/s');
% Calculate mean and standard deviation, and display mean and ± two standard deviations with horizontal lines Velocity_Mean = mean(Velocity_Instantaneous); hold on; plot([t(1) t(end)], [1 1]*Velocity_Mean,'--k') Velocity_STD = std(Velocity_Instantaneous); plot([t(1) t(end)], [1 1]*(Velocity_Mean-Velocity_STD*2),':m','linewidth',1.2) plot([t(1) t(end)], [1 1]*(Velocity_Mean+Velocity_STD*2),':m','linewidth',1.2) Turbulence_Intensity = Velocity_STD/Velocity_Mean;
% Calculate the Skewness and Kurtosis of the data %'skewness' and kurtosis requires Statistics and Machine Learning Toolbox.% Velocity_SK = skewness(Velocity_Instantaneous)% Velocity_KU = kurtosis(Velocity_Instantaneous)
% Investigate the effect of sampling duration on the mean j=1; block_length=100:100:4000; for i=block_length velocity_std_diff_length(j)=std(Velocity_Instantaneous(1:i)); j=j+1; end figure, plot(t(block_length),velocity_std_diff_length,'o-');grid on; axis([0 1 0.8 1.2]) xlabel('T, s');ylabel('\sigma(T), m/s');%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Calculate frequency spectrumPower_Spectrum_Mean=zeros(Sample_Number_Per_Block,1);for Block_Selected=1:Number_Blocks Velocity_Instantaneous=Velocity_Loaded(1:Sample_Number_Per_Block,Block_Selected); Velocity_Mean = mean(Velocity_Instantaneous); Velocity_deMean=Velocity_Instantaneous-Velocity_Mean;
Velocity_FFT=fft(Velocity_deMean); Power_Spectrum=Velocity_FFT.*conj(Velocity_FFT)/(Sample_Number_Per_Block*Sample_Frequency); Power_Spectrum_Mean=Power_Spectrum_Mean+Power_Spectrum;endPower_Spectrum_Mean=Power_Spectrum_Mean/Number_Blocks;
% Check if the integral of the spectrum equals the variance of the original data, the value of the following equation should be approximately equal to 1sum(Power_Spectrum)*df/var(Velocity_Instantaneous)
fft_n=Sample_Number_Per_Block/2+1;PSD_OneSide=Power_Spectrum(1:fft_n)*2;PSD_OneSide_Mean=Power_Spectrum_Mean(1:fft_n)*2;
figure, loglog(f,Power_Spectrum,'b')hold on;loglog(f,Power_Spectrum_Mean,'r');xlabel('f, Hz');ylabel('E_{uu}, (m/s)^2/Hz');title('Two-sided frequency spectrum')
figure, loglog(f(1:fft_n),PSD_OneSide);grid on;hold on;loglog(f(1:fft_n),PSD_OneSide_Mean,'r');f_sp=[100 1000]; %Frequency, Hzsp=f_sp.^(-5/3)*100; %loglog(f_sp,sp,'-m');xlabel('f, Hz');ylabel('E_{uu}, (m/s)^2/Hz');legend('Spectrum using 1 block','Spectrum averaged from 180 blocks','-5/3','location','SouthWest')title('One-sided frequency spectrum')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calculate wavenumber spectrumk_x=2*pi*f/Velocity_Mean;Wave_PSD_OneSide_Mean=PSD_OneSide_Mean*Velocity_Mean/(2*pi);figure, loglog(k_x(1:fft_n),Wave_PSD_OneSide_Mean.*k_x(1:fft_n)','r');grid on;xlabel('k_x, rad');ylabel('\Phi_{uu}*k_x, (m/s)^2');title('Pre-multiplied wave spectrum')k_t=[100 2000];P_t=k_t.^(-5/3)*10;hold on; loglog(k_t,P_t.*k_t,'-k');
%%%% Calculate minimum scale, unitmdk=k_x(2)-k_x(1);Kine_Viscosity=1.51e-5; % Kinematic viscosity of air at 20 Celsius, m/s^2Strain_Rate_ISO=sum(PSD_OneSide_Mean.*(k_x(1:fft_n)'.^2))*dk; Epsilon_ISO=15*Kine_Viscosity*Strain_Rate_ISO; % dissipation rateKolm_Scale=(Kine_Viscosity^3/Epsilon_ISO)^.25
% Calculate Taylor scale, unitmVelocity_Variance=mean(var(Velocity_Loaded));Taylor_Scale=sqrt(Velocity_Variance/Strain_Rate_ISO)
% Calculate Reynolds number based on Taylor scaleRe_Lambda=sqrt(Velocity_Variance)*Taylor_Scale/Kine_Viscosity%C_e=Epsilon_ISO*.0183/Velocity_Variance.^1.5
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Create time difference arraytau=[-Sample_Number_Per_Block+1:1:Sample_Number_Per_Block-1]/Sample_Frequency;
% Calculate autocorrelation coefficientAutoCorrelation_Mean=zeros(Sample_Number_Per_Block*2-1,1);for Block_Selected=1:Number_Blocks Velocity_Instantaneous=Velocity_Loaded(1:Sample_Number_Per_Block,Block_Selected); Velocity_Mean = mean(Velocity_Instantaneous); Velocity_deMean=Velocity_Instantaneous-Velocity_Mean;
AutoCorrelation = xcorr(Velocity_deMean); AutoCorrelation = AutoCorrelation/max(AutoCorrelation); AutoCorrelation_Mean = AutoCorrelation_Mean + AutoCorrelation;endAutoCorrelation_Mean = AutoCorrelation_Mean/Number_Blocks;
% Plot the complete autocorrelation functionfigure, plot(tau,AutoCorrelation_Mean);grid on;xlabel('\tau, s');ylabel('\rho_{uu}(\tau)');
% Plot the one-sided autocorrelation function and zoom in locally.% Determine the time difference \tau when the function crosses the x-axis and mark it with a red circle.d_tau=1/Sample_Frequency;AutoCorrelation_Mean_Right_Half=AutoCorrelation_Mean(Sample_Number_Per_Block:end);tau_Right_Half=tau(Sample_Number_Per_Block:end);[temp_a,temp_b]=findpeaks(-abs(AutoCorrelation_Mean_Right_Half));figure, plot(tau_Right_Half,AutoCorrelation_Mean_Right_Half);hold on;grid on;plot(tau_Right_Half(temp_b(1)),AutoCorrelation_Mean_Right_Half(temp_b(1)),'or');hold on; plot([0 .1],[0 0],'k-');axis([0 .05 -.2 1])
% Calculate integral time scale and length scale, units are s and mIntegral_Time_Scale=sum(AutoCorrelation_Mean_Right_Half(1:temp_b(1)))*d_tauIntegral_Length_Scale=Integral_Time_Scale*Velocity_Mean
% Calculate the ratio of sampling duration to 2 times the integral scale to determine the number of independent samplesNs=T/(2*Integral_Time_Scale);
% Use the number of independent samples to estimate the measurement error of mean and standard deviationepsilon_mean=std(Velocity_Instantaneous)/mean(Velocity_Instantaneous)/sqrt(Ns)epsilon_std=1/sqrt(Ns)
