PowerSpecutrum.m

% Compute power spectrum

% Generated by MATLAB(R) 9.11 and Signal Processing Toolbox 8.7.
% Generated on: 10-Sep-2022 16:24:48

%%
% Power spectral analysis
% Add path of load_open_ephys_binary to Matlab first.
addpath = ('/Users/zhangjinming/Desktop/MatlabCode')
D = load_open_ephys_binary(['/Users/zhangjinming/Desktop/实验结果/博士课题/多通道电生理' ...
    '/JZL195/JZL195/24h after i.p./2#-32Ch_2022-08-31_13-51-14_24h after JZL195/Record Node 113/experiment1/recording1/structure.oebin'],'continuous',1);
for r = 17:32 %Load raw data per channel (16-channel)
    eval(['Y' num2str(r-16) '=D.Data(r,:)']) 
end 

%批量降低采样评率30kHz==500Hz 
for e=1:16
 eval(['Y_' num2str(e)  '=resample(Y' num2str(e) ',500, 30000);']); 
end

% Parameters
timeLimits = [-28.43998 653.5401]; % seconds
frequencyLimits = [0.5 150]; % Hz

%%
% Index into signal time region of interest
Y_1_ROI = Y_1(:);
sampleRate = 500; % Hz
startTime = 0; % seconds
minIdx = ceil(max((timeLimits(1)-startTime)*sampleRate,0))+1;
maxIdx = floor(min((timeLimits(2)-startTime)*sampleRate,length(Y_1_ROI)-1))+1;
Y_1_ROI = normalize(Y_1_ROI(minIdx:maxIdx));


% Compute spectral estimate
% Run the function call below without output arguments to plot the results
[PY_1_ROI, FY_1_ROI] = pspectrum(Y_1_ROI,sampleRate, ...
    'FrequencyLimits',frequencyLimits);

%%
% Index into signal time region of interest
Y_10_ROI = Y_10(:);
sampleRate = 500; % Hz
startTime = 0; % seconds
minIdx = ceil(max((timeLimits(1)-startTime)*sampleRate,0))+1;
maxIdx = floor(min((timeLimits(2)-startTime)*sampleRate,length(Y_10_ROI)-1))+1;
Y_10_ROI = normalize(Y_10_ROI(minIdx:maxIdx));

% Compute spectral estimate
% Run the function call below without output arguments to plot the results
[PY_10_ROI, FY_10_ROI] = pspectrum(Y_10_ROI,sampleRate, ...
    'FrequencyLimits',frequencyLimits);

%%
% Index into signal time region of interest
Y_11_ROI = Y_11(:);
sampleRate = 500; % Hz
startTime = 0; % seconds
minIdx = ceil(max((timeLimits(1)-startTime)*sampleRate,0))+1;
maxIdx = floor(min((timeLimits(2)-startTime)*sampleRate,length(Y_11_ROI)-1))+1;
Y_11_ROI = normalize(Y_11_ROI(minIdx:maxIdx));

% Compute spectral estimate
% Run the function call below without output arguments to plot the results
[PY_11_ROI, FY_11_ROI] = pspectrum(Y_11_ROI,sampleRate, ...
    'FrequencyLimits',frequencyLimits);

%%
% Index into signal time region of interest
Y_12_ROI = Y_12(:);
sampleRate = 500; % Hz
startTime = 0; % seconds
minIdx = ceil(max((timeLimits(1)-startTime)*sampleRate,0))+1;
maxIdx = floor(min((timeLimits(2)-startTime)*sampleRate,length(Y_12_ROI)-1))+1;
Y_12_ROI = normalize(Y_12_ROI(minIdx:maxIdx));

% Compute spectral estimate
% Run the function call below without output arguments to plot the results
[PY_12_ROI, FY_12_ROI] = pspectrum(Y_12_ROI,sampleRate, ...
    'FrequencyLimits',frequencyLimits);

%%
% Index into signal time region of interest
Y_13_ROI = Y_13(:);
sampleRate = 500; % Hz
startTime = 0; % seconds
minIdx = ceil(max((timeLimits(1)-startTime)*sampleRate,0))+1;
maxIdx = floor(min((timeLimits(2)-startTime)*sampleRate,length(Y_13_ROI)-1))+1;
Y_13_ROI = normalize(Y_13_ROI(minIdx:maxIdx));

% Compute spectral estimate
% Run the function call below without output arguments to plot the results
[PY_13_ROI, FY_13_ROI] = pspectrum(Y_13_ROI,sampleRate, ...
    'FrequencyLimits',frequencyLimits);

%%
% Index into signal time region of interest
Y_14_ROI = Y_14(:);
sampleRate = 500; % Hz
startTime = 0; % seconds
minIdx = ceil(max((timeLimits(1)-startTime)*sampleRate,0))+1;
maxIdx = floor(min((timeLimits(2)-startTime)*sampleRate,length(Y_14_ROI)-1))+1;
Y_14_ROI = normalize(Y_14_ROI(minIdx:maxIdx));

% Compute spectral estimate
% Run the function call below without output arguments to plot the results
[PY_14_ROI, FY_14_ROI] = pspectrum(Y_14_ROI,sampleRate, ...
    'FrequencyLimits',frequencyLimits);

%%
% Index into signal time region of interest
Y_15_ROI = Y_15(:);
sampleRate = 500; % Hz
startTime = 0; % seconds
minIdx = ceil(max((timeLimits(1)-startTime)*sampleRate,0))+1;
maxIdx = floor(min((timeLimits(2)-startTime)*sampleRate,length(Y_15_ROI)-1))+1;
Y_15_ROI = normalize(Y_15_ROI(minIdx:maxIdx));

% Compute spectral estimate
% Run the function call below without output arguments to plot the results
[PY_15_ROI, FY_15_ROI] = pspectrum(Y_15_ROI,sampleRate, ...
    'FrequencyLimits',frequencyLimits);

%%
% Index into signal time region of interest
Y_16_ROI = Y_16(:);
sampleRate = 500; % Hz
startTime = 0; % seconds
minIdx = ceil(max((timeLimits(1)-startTime)*sampleRate,0))+1;
maxIdx = floor(min((timeLimits(2)-startTime)*sampleRate,length(Y_16_ROI)-1))+1;
Y_16_ROI = normalize(Y_16_ROI(minIdx:maxIdx));

% Compute spectral estimate
% Run the function call below without output arguments to plot the results
[PY_16_ROI, FY_16_ROI] = pspectrum(Y_16_ROI,sampleRate, ...
    'FrequencyLimits',frequencyLimits);

%%
% Index into signal time region of interest
Y_2_ROI = Y_2(:);
sampleRate = 500; % Hz
startTime = 0; % seconds
minIdx = ceil(max((timeLimits(1)-startTime)*sampleRate,0))+1;
maxIdx = floor(min((timeLimits(2)-startTime)*sampleRate,length(Y_2_ROI)-1))+1;
Y_2_ROI = normalize(Y_2_ROI(minIdx:maxIdx));

% Compute spectral estimate
% Run the function call below without output arguments to plot the results
[PY_2_ROI, FY_2_ROI] = pspectrum(Y_2_ROI,sampleRate, ...
    'FrequencyLimits',frequencyLimits);

%%
% Index into signal time region of interest
Y_3_ROI = Y_3(:);
sampleRate = 500; % Hz
startTime = 0; % seconds
minIdx = ceil(max((timeLimits(1)-startTime)*sampleRate,0))+1;
maxIdx = floor(min((timeLimits(2)-startTime)*sampleRate,length(Y_3_ROI)-1))+1;
Y_3_ROI = normalize(Y_3_ROI(minIdx:maxIdx));

% Compute spectral estimate
% Run the function call below without output arguments to plot the results
[PY_3_ROI, FY_3_ROI] = pspectrum(Y_3_ROI,sampleRate, ...
    'FrequencyLimits',frequencyLimits);

%%
% Index into signal time region of interest
Y_4_ROI = Y_4(:);
sampleRate = 500; % Hz
startTime = 0; % seconds
minIdx = ceil(max((timeLimits(1)-startTime)*sampleRate,0))+1;
maxIdx = floor(min((timeLimits(2)-startTime)*sampleRate,length(Y_4_ROI)-1))+1;
Y_4_ROI = normalize(Y_4_ROI(minIdx:maxIdx));

% Compute spectral estimate
% Run the function call below without output arguments to plot the results
[PY_4_ROI, FY_4_ROI] = pspectrum(Y_4_ROI,sampleRate, ...
    'FrequencyLimits',frequencyLimits);

%%
% Index into signal time region of interest
Y_5_ROI = Y_5(:);
sampleRate = 500; % Hz
startTime = 0; % seconds
minIdx = ceil(max((timeLimits(1)-startTime)*sampleRate,0))+1;
maxIdx = floor(min((timeLimits(2)-startTime)*sampleRate,length(Y_5_ROI)-1))+1;
Y_5_ROI = normalize(Y_5_ROI(minIdx:maxIdx));

% Compute spectral estimate
% Run the function call below without output arguments to plot the results
[PY_5_ROI, FY_5_ROI] = pspectrum(Y_5_ROI,sampleRate, ...
    'FrequencyLimits',frequencyLimits);

%%
% Index into signal time region of interest
Y_6_ROI = Y_6(:);
sampleRate = 500; % Hz
startTime = 0; % seconds
minIdx = ceil(max((timeLimits(1)-startTime)*sampleRate,0))+1;
maxIdx = floor(min((timeLimits(2)-startTime)*sampleRate,length(Y_6_ROI)-1))+1;
Y_6_ROI = normalize(Y_6_ROI(minIdx:maxIdx));

% Compute spectral estimate
% Run the function call below without output arguments to plot the results
[PY_6_ROI, FY_6_ROI] = pspectrum(Y_6_ROI,sampleRate, ...
    'FrequencyLimits',frequencyLimits);

%%
% Index into signal time region of interest
Y_7_ROI = Y_7(:);
sampleRate = 500; % Hz
startTime = 0; % seconds
minIdx = ceil(max((timeLimits(1)-startTime)*sampleRate,0))+1;
maxIdx = floor(min((timeLimits(2)-startTime)*sampleRate,length(Y_7_ROI)-1))+1;
Y_7_ROI = normalize(Y_7_ROI(minIdx:maxIdx));

% Compute spectral estimate
% Run the function call below without output arguments to plot the results
[PY_7_ROI, FY_7_ROI] = pspectrum(Y_7_ROI,sampleRate, ...
    'FrequencyLimits',frequencyLimits);

%%
% Index into signal time region of interest
Y_8_ROI = Y_8(:);
sampleRate = 500; % Hz
startTime = 0; % seconds
minIdx = ceil(max((timeLimits(1)-startTime)*sampleRate,0))+1;
maxIdx = floor(min((timeLimits(2)-startTime)*sampleRate,length(Y_8_ROI)-1))+1;
Y_8_ROI = normalize(Y_8_ROI(minIdx:maxIdx));

% Compute spectral estimate
% Run the function call below without output arguments to plot the results
[PY_8_ROI, FY_8_ROI] = pspectrum(Y_8_ROI,sampleRate, ...
    'FrequencyLimits',frequencyLimits);

%%
% Index into signal time region of interest
Y_9_ROI = Y_9(:);
sampleRate = 500; % Hz
startTime = 0; % seconds
minIdx = ceil(max((timeLimits(1)-startTime)*sampleRate,0))+1;
maxIdx = floor(min((timeLimits(2)-startTime)*sampleRate,length(Y_9_ROI)-1))+1;
Y_9_ROI = normalize(Y_9_ROI(minIdx:maxIdx));

% Compute spectral estimate
% Run the function call below without output arguments to plot the results
[PY_9_ROI, FY_9_ROI] = pspectrum(Y_9_ROI,sampleRate, ...
    'FrequencyLimits',frequencyLimits);



PY(:,1) = PY_1_ROI(:,1);
PY(:,2) = PY_2_ROI(:,1);
PY(:,3) = PY_3_ROI(:,1);
PY(:,4) = PY_4_ROI(:,1);
PY(:,5) = PY_5_ROI(:,1);
PY(:,6) = PY_6_ROI(:,1);
PY(:,7) = PY_7_ROI(:,1);
PY(:,8) = PY_8_ROI(:,1);
PY(:,9) = PY_9_ROI(:,1);
PY(:,10) = PY_10_ROI(:,1);
PY(:,11) = PY_11_ROI(:,1);
PY(:,12) = PY_12_ROI(:,1);
PY(:,13) = PY_13_ROI(:,1);
PY(:,14) = PY_14_ROI(:,1);
PY(:,15) = PY_15_ROI(:,1);
PY(:,16) = PY_16_ROI(:,1);