Emotion-Detection/STFT.m at master · xinzhu-cai/Emotion-Detection · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
function [ psd, de ] = STFT(data,stft_para)
%input: data [n*m]          n electrodes, m time points
%       stft_para.stftn     frequency domain sampling rate
%       stft_para.fStart    start frequency of each frequency band
%       stft_para.fEnd      end frequency of each frequency band
%       stft_para.window    window length of each sample point(seconds)
%       stft_para.f         original frequency
%output:psd,DE [n*l*k]        n electrodes, l windows, k frequency bands

    %initialize the parameters
    STFTN=stft_para.stftn;
    fStart=stft_para.fStart;
    fEnd=stft_para.fEnd;
    fs=stft_para.fs;
    window=stft_para.window;

    WindowPoints=fs*window;

    %%Transform the frequency from original one to sampled
    fStartNum=zeros(1,length(fStart));
    fEndNum=zeros(1,length(fEnd));
    for i=1:length(fStart)
        fStartNum(1,i)=fix(fStart(i)/fs*STFTN);
        fEndNum(1,i)=fix(fEnd(i)/fs*STFTN)-1;
    end

    [n m]=size(data);
    % Number of windows
    l=fix(m/WindowPoints);
    psd=zeros(n,l,length(fStart));
    de = zeros(n, l, length(fStart));
    %get Hanning window function
    Hlength=window*fs;
    Hwindow=hanning(Hlength);
    %New window length
    WindowPoints=fs*window;
    for i=1:l
        dataNow=data(:,WindowPoints*(i-1)+1:WindowPoints*i);
        for j=1:n
            %For each electrode
            temp=dataNow(j,:);
            % x[n] * w[n-m]
            Hdata=temp.*Hwindow';
            % STFTN is number of points
            FFTdata=fft(Hdata,STFTN);
            % |X(m,w)|
            magFFTdata=abs(FFTdata(1:1:STFTN/2));
            for p=1:length(fStart)
                E=0;
                E_log = 0;
                for p0=fStartNum(p):fEndNum(p)
                    E=E+magFFTdata(p0)*magFFTdata(p0);
                %    E_log = E_log + log2(magFFTdata(p0)*magFFTdata(p0)+1);
                end
                E=E/(fEndNum(p)-fStartNum(p)+1);
                psd(j,i,p)= E;
                de(j,i,p) = log2(100*E);
                %de(j,i,p)=log2((1+E)^4);
            end
        end
    end
end