Homec4science
Diffusion PP routine (master) Browse

R11173/Sxx_Function_Local.m
f2cd617ca124master

Last Change
Raw File

Sxx_Function_Local.m
View Options

%% Arthur Favrel + Christian Landry (2015, Hyperbole project)
function [Frequency, Sxx] = Sxx_Function_Local(Time, Signal, Method, Overlap, Window_Size)
NFFT = 2^nextpow2(Window_Size);
%NFFT=Window_Size;
L = length(Time); % Signal Size
dt = Time(2) - Time(1); % resolution dt
if strcmp(Method, 'Hann')
Windowing = hann(NFFT);
Constante = 2;
elseif strcmp(Method, 'Hamming')
Windowing = hamming(NFFT);
Constante = 1/mean(Windowing); % the correction factor is eqaul to the inverse of the mean of the weight factor of the window. = env 1.85 for Hamming
elseif strcmp(Method, 'Rectangle')
Windowing = rectwin(NFFT);
Constante = 1;
end
%% Subdivision of the signal signal
overlapp = ceil(Overlap*NFFT);
beg = []; % Vector to store first indice of each sub-signal
% First sub-division
beg = [beg 1];
% Next sub-division
j = 1;
while beg(j) + 2*NFFT - overlapp - 1 < L
beg = [beg beg(j) + NFFT - overlapp];
j = j + 1;
end
Nd = length(beg); % Maximum number of sub-signal with overlapping
for j = 1:Nd
Signal_Sub_1(:, j) = Signal(beg(j):beg(j)+NFFT-1);
end
%% Statistical value computation - on each sub-signal
% Mean value
mean1 = [];
for j = 1:Nd
mean1(j) = mean(Signal_Sub_1(:, j));
end
%% Centered Sub-signals
for j = 1:Nd
Signal_Sub_centered_1(:, j) = Signal_Sub_1(:, j) - mean1(j);
end
%% Double sided spectra computation on each sub-signal
% Double sided spectra - initial signal (The length signal are taken 2^n for this case)
frequency_double = [0:(NFFT-1)]*(1/(dt*NFFT));
% With windowing
for j = 1:Nd
fft_signal_shared(:, j) = abs(Constante*fft(Windowing.*Signal_Sub_centered_1(:,j)));
end
% Average all windows
fxx_av = zeros(NFFT, 1);
for j=1:Nd
fxx_av = fxx_av + fft_signal_shared(:,j);
end
fxx_av = fxx_av./Nd;
% Standard procedure for matlab fft (See doc fft)
P2=abs(fxx_av/NFFT);
P1=P2(1:NFFT/2+1);
P1(2:end-1) = 2*P1(2:end-1);
Sxx = P1;
Frequency = frequency_double(1:NFFT/2+1); %See definition of one-sided spectra
c4science · Help