I want to do loop operation of these files but fail to get the desired result. Can anyone help me to figure out the problem?

After the operation, I wish to get the values as in the table:

The script that I work on, can't store/save the data as in the table format.

p = dir('f*.txt');
nFiles = numel(p);
counter=0;
for kk = 1:nFiles
        counter= kk+1;
        fileID1 = fopen(p(kk).name,'rt');
        C_data = textscan(fileID1,'%f');
        f1 = cell2mat(C_data);  % A = dlmread(p(k).name);  
        %A_cat{kk} = A; disp(size(A_cat));
              fclose(fileID1);
            val = cell(1, nFiles-1); % Memory Pre-allocation
            for i =  counter: nFiles
                fileID2 = fopen(p(i).name,'rt');
                C_data = textscan(fileID2,'%f');
                f2 = cell2mat(C_data);  % A = dlmread(p(k).name);
                fclose(fileID2);
                [Pyy,freq]=cpsd(f2,f2,hanning(512),[],512,500);
                [Pxy,freq]=cpsd(f1,f2,hanning(512),[],512,500);
                [Pxx,freq]=cpsd(f1,f1,hanning(512),[],512,500);
                coh= (Pxy)./sqrt(Pxx.*Pyy);
                val=real(coh(42))
                %val=horzcat(val) 
            end
           val(kk)=horzcat(val) 
end

Can anyone help to fix the issues? Thank you in advance.

Answer 1

[Mathieu NOE]

hello

I tweaked a bit your code , if fact it was simpler in my opinion to do the full 4 x 4 loops instead of doing the specific counter condition

now there is one point that surprises me, its the diagonal of your matrice . the autocorrelation of each element with itself should be 1 and not zero (and it's what I get) - so I wonder about this point. of course we can change the code to force the diagonal to be zero but that sounds strange to me

p = dir('f*.txt');
nFiles = numel(p);
counter=0;
for kk = 1:nFiles
        %counter= kk+1
        fileID1 = fopen(p(kk).name,'rt');
        C_data = textscan(fileID1,'%f');
        f1 = cell2mat(C_data);  % A = dlmread(p(k).name);  
        %A_cat{kk} = A; disp(size(A_cat));
              fclose(fileID1);
            %val = cell(1, nFiles-1); % Memory Pre-allocation
            for i =  1: nFiles
                fileID2 = fopen(p(i).name,'rt');
                C_data = textscan(fileID2,'%f');
                f2 = cell2mat(C_data);  % A = dlmread(p(k).name);
                fclose(fileID2);
                [Pyy,freq]=cpsd(f2,f2,hanning(512),[],512,500);
                [Pxy,freq]=cpsd(f1,f2,hanning(512),[],512,500);
                [Pxx,freq]=cpsd(f1,f1,hanning(512),[],512,500);
                coh= (Pxy)./sqrt(Pxx.*Pyy);
                val(kk,i)=real(coh(42)) ;
            end
end

% gives : 

val =

    1.0000    0.6725   -0.2602    0.3779
    0.6725    1.0000    0.1127    0.2854
   -0.2602    0.1127    1.0000    0.0763
    0.3779    0.2854    0.0763    1.0000
    

[SA]

@Mathieu NOE, Thank you very much for explaining the issue. Obviously, you are right (diagonal should be 1, not 0) and it works fine.

When I do memory pre-allocation

val = cell(nFiles); % Memory Pre-allocation 
......
    val{kk,i}=real(coh(42)) ;    

script executes but return an empty cell array. How to check the cell array 'val'? Thanks.

[Mathieu NOE]

hello again

the size in the preallocation is not corret

also do it before the for loop

p = dir('f*.txt');
nFiles = numel(p);
counter=0;            
val = cell(nFiles, nFiles); % Memory Pre-allocation HERE

for kk = 1:nFiles
        %counter= kk+1
        fileID1 = fopen(p(kk).name,'rt');
        C_data = textscan(fileID1,'%f');
        f1 = cell2mat(C_data);  
              fclose(fileID1);
            for i =  1: nFiles
                fileID2 = fopen(p(i).name,'rt');
                C_data = textscan(fileID2,'%f');
                f2 = cell2mat(C_data);
                fclose(fileID2);
                [Pyy,freq]=cpsd(f2,f2,hanning(512),[],512,500);
                [Pxy,freq]=cpsd(f1,f2,hanning(512),[],512,500);
                [Pxx,freq]=cpsd(f1,f1,hanning(512),[],512,500);
                coh= (Pxy)./sqrt(Pxx.*Pyy);
                val{kk,i}=real(coh(42)) ;
            end
end

[SA]

Hi,

Again, I come here. I tried with this code but the process is too slow to execute. Do you have any idea to speed up the process? In the 'val' (symmetric) matrix if I calculate only upper/lower value from the diagonal then I may save around 50% of time as

val{i,j}=val{j,i}

another point, as I don't use 'freq' so I've replaced 'freq' by '~'

[Pyy,~]=cpsd(f2,f2,hanning(512),[],512,500)

N.B.: I've many files and size of each files 50M bytes. Thank you.

[Mathieu NOE]

hello

so you are doing a cross correlation matrix based on time series , that's what I understand today

why are the data files so big ? are you sure that you need all this amount of data to process, for example maybe your files at sampled at too high sampling freq and this could be reduced to speed up the process (using decimate); have you conducted a king of frequency analysis to define in which frequency range the data must be ?

second, why do you pick only the 42th value of the coherence ? (coh(42))

[SA]

Thanks for your question. File size is big as I consider longer time (10,000+ seconds). I want to check the Real Coherence value at a specific frequency(say @40Hz =(coh(42))). Here, the sampling frequency is chosen fs=500 and 'noverlap' left blank which matlab set 50% (default).

so, I tried the following scripts.

clc; close all;
p = dir('f*.txt');
nFiles = numel(p);
counter=0;
val = zeros(nFiles); % Memory Pre-allocation
for kk = 1:nFiles
        %counter= kk+1
        fileID1 = fopen(p(kk).name,'rt');
        C_data = textscan(fileID1,'%f');
        f1 = cell2mat(C_data);  % A = dlmread(p(k).name);  
        fclose(fileID1);
            
            for i =  (kk+1): nFiles
                fileID2 = fopen(p(i).name,'rt');
                C_data = textscan(fileID2,'%f');
                f2 = cell2mat(C_data);  % A = dlmread(p(k).name);
                fclose(fileID2);
                [Pyy,~]=cpsd(f2,f2,hanning(512),[],512,500);
                [Pxy,~]=cpsd(f1,f2,hanning(512),[],512,500);
                [Pxx,~]=cpsd(f1,f1,hanning(512),[],512,500);
                coh= (Pxy)./sqrt(Pxx.*Pyy);
                realCoh=real(coh(42));
                val(kk,i)=val(i,kk)+realCoh;
            end
end
final_val=val+val'+eye(nFiles);

If I calculate the half of the symmetric matrix and adding transpose of it, it saves 60% of time. Thanks for your time.