Calculating trembling frequency

calculating_trembling_frequency.m

%% (0) Set option
vid_save_option = false; % <- change this to true (to gen video output)

%% (1) Read frames from input video
vid_in = VideoReader( 'video_input.mp4' );
frames = [];
nFrames = vid_in.Duration * vid_in.FrameRate;
for frameIdx = 1:nFrames
    frames(:,:,:,frameIdx) = vid_in.readFrame();
end

%% (2) Prepare output video object
if vid_save_option
    try, vid_out.close(); end
    vid_out = VideoWriter('video_output', 'MPEG-4');
    vid_out.FrameRate = round(vid_in.FrameRate);
    vid_out.open();
end

%% (3) Set sample coordinate (target to detect)
coordinate = [];
win_x = [210:340];
win_y = [200:350];
shoe = (frames(win_x,win_y,:,1));
target_to_detect = rgb2gray(uint8(shoe));

%% (4) Open figure handle
figure_to_draw = figure(2); clf;
set(gcf, 'Color', 'w', 'Position', [0 200 900 500]);
fontsize = 10;

%% (5) Main loop
for frameIdx = 1:nFrames
    
    % Draw raw frame
    subplot(2,2,1);
    img = frames(:,:,:,frameIdx);
    imagesc( uint8( img ))
    tt = round(1000*frameIdx/30)/1000;
    title(['Raw data, t = ' sprintf( '%.3f' , tt) 's']);
    set(gca, 'FontSize', fontsize, 'LineWidth', 2, 'Box', 'off' );
    grid on
    xlabel('X (unit: pixel)'); ylabel('Y (unit: pixel)');
    
    % Draw template to detect
    if frameIdx==1
        subplot(2,2,2);
        imshow( uint8( shoe ));
        title([ 'Template: ' sprintf( 'x=[%d~%d]', win_x([1,end])) ', ' sprintf( 'y=[%d~%d]', win_y([1,end])) ]);
        set(gca, 'FontSize', fontsize, 'LineWidth', 2, 'Box', 'off' );
    end
    
    % Detect target
    grayframe = rgb2gray(uint8(img));%rgb2gray(imread('peppers.png'));
    c = normxcorr2(target_to_detect,grayframe);
    [ypeak, xpeak] = find(c==max(c(:)));
    yoffSet = ypeak-size(target_to_detect,1);
    xoffSet = xpeak-size(target_to_detect,2);
    
    % Save time series data
    coordinate(1,frameIdx) = xoffSet;
    coordinate(2,frameIdx) = yoffSet;
    
    % Show location of detected target (demonstration)
    subplot(2,2,3); hold off;
    imagesc(uint8(img));
    rect = imrect(gca, ...
        [xoffSet+1, yoffSet+1, ...
        size(target_to_detect,2), size(target_to_detect,1)] );
    rect.setColor([0 1 0]);
    xlabel('X (unit: pixel)'); ylabel('Y (unit: pixel)');
    title('Detected Shoe');
    set(gca, 'FontSize', fontsize, 'LineWidth', 2, 'Box', 'off' );
    hold on
    plot( [0 1]*xoffSet, [1 1]*yoffSet,  'g-' )
    text( .5*xoffSet, yoffSet-30, sprintf('Y = %3d', yoffSet), 'Color', 'g')
    
    %% Plot time domain data
    subplot(2,4,7); hold off;
    time_vector = 1000*[1:frameIdx]/30;
    plot( time_vector, coordinate(2,:), 'k' );
    text( time_vector(end)+10, coordinate(2,end), ...
        sprintf('Y = %3d', yoffSet), 'Color', 'g' );
    
    title('Time domain')
    xlabel('Time (msec)')
    ylabel('Y position (pixel)');
    set(gca, 'YDir', 'reverse')
    xlim([0 4000]);
    %     ylim([180 250]);
    ylim([1 404]);
    set(gca, 'FontSize', fontsize, 'LineWidth', 2, 'Box', 'off' );
    
    % Plot freq-domain data
    if frameIdx > 10
        subplot(2,4,8); hold off;
        [X,freq]=positiveFFT( coordinate(2,:), 30 );
        power = abs(X).^2;
        plot( freq, power, 'ko-' );
        xlim([1 10]);
        ylim([0 100])
        xlabel('Freq (Hz)')
        ylabel('Power (pixel^2/Hz)');
        title(['Frequency domain'])
        set(gca, 'FontSize', fontsize, 'LineWidth', 2, 'Box', 'off' );
        set(gca, 'XTick', 0:1:100, 'YTick', 0:20:100)
        
        hold on;
        % Peak detection
        freq_inspection_range = [2, 10]; % Hz
        freq_idx = hb_findIdx( freq_inspection_range, freq );
        [peak_value,max_idx] = max( power(freq_idx) );
        peak_hz = freq( freq_idx(max_idx) );
        
        plot( peak_hz, peak_value, 'r*' )
        text( peak_hz+1, peak_value, sprintf( 'Peak: %.1f Hz', peak_hz) );
    end
    
%     % Add title
%     if frameIdx==1
%         try
%             spl=suptitle(['Theta Rhythmicity of Human Shoe']);
%             spl.Position = [0.5000 -0.022 0];
%             spl.FontSize = 14;
%         end
%     end
    
    % Draw 
    drawnow;
    if vid_save_option
        f=getframe(figure_to_draw);
        vid_out.writeVideo(f.cdata);
    end
    
end
try, vid_out.close(); end

%% Subfunctions
function [idx,short]=hb_findIdx(range,fullData)
idx=max(find(fullData<range(1)))+1:max(find(fullData<range(2)));
end
function [X,freq]=positiveFFT(x,Fs)
N=(length(x));  %get the number of points
k=0:N-1;        %create a vector from 0 to N-1
T=N/Fs;         %get the frequency interval
freq=k/T;       %create the frequency range
cutOff = ceil(N/2);
freq = freq(1:cutOff);
X=fft(x)/N*2; % normalize the data
X = X(1:cutOff); % Single trial FFT
end