Animated magnitude spectrum (windowed fft)

Hi, I need to show the Magnitude spectrum made up of a window of 128 samples for the loaded signals. These two-magnitude spectra need to be animated (plotted) w.r.t. the time-domain plots on the LHS of the figure. Till now I have manged to extract the magnitude spectrum for the whole signal and then plot after the animation, any ideas how I can alter the code so that it doses the above mentioned? Thanks!

clear all
close all
clc

load('ch1_first_5s.mat');
load('ch2_first_5s.mat');

fs = 128;
n1 = 0:1:((5 * fs) - 1);
n2 = 1:1:(fs - 1);

for i = 1:1:511
    
    subplot(2, 2, 1);
    plot((i + n2), ch1_first_5s(i + n2), 'b');
    title('Channel 1 - Time Domain EEG plot');
    ylim([-200 200]); ylabel('Amplitude');
    xlim([i (i + 128)]); xlabel('Discrete Time n (Sample Number)');
    grid on;

    subplot(2, 2, 3);
    plot((i + n2), ch2_first_5s(i + n2), 'b');
    title('Channel 2 - Time Domain EEG plot');
    ylim([-200 200]); ylabel('Amplitude');
    xlim([i (i + 128)]); xlabel('Discrete Time n (Sample Number)');
    grid on;
        
    pause(0.05);

end

N = length(n1);
N1 = (2)^nextpow2(N);
X_k1 = fft(ch1_first_5s, N1); X_k1 = X_k1(1:(N1 / 2));
X_k2 = fft(ch2_first_5s, N1); X_k2 = X_k2(1:(N1 / 2));
Mag1 = abs(X_k1);
Mag2 = abs(X_k2);
f = fs * (0:(N1 / 2) - 1) / N1; 

subplot(2, 2, 2);
plot(f, mag2db(Mag1 / N1), 'b');
ylim([-60 20]); ylabel('Magnitude (dB)');
xlim([0 f(end)]); xlabel('Frequency (Hz)');
title('Channel 1 - Freq Domain EEG plot');
grid on;

subplot(2, 2, 4);
plot(f, mag2db(Mag2 / N1), 'b'); 
ylim([-60 20]); ylabel('Magnitude (dB)');
xlim([0 f(end)]); xlabel('Frequency (Hz)');
title('Channel 2 - Freq Domain EEG plot');
grid on;

Answer 1

[Mathieu NOE]

hello Joshua

this is my suggestion / code below.

It shows how to update the contents of a plot without (re)creating the entire plot structure at each time iteration (which can be pretty slow for complex figures)

so the plot is created only once and then we update some data and axes properties / values

I tested my code also on dummy sine waves just to make sure the correct data is displayed in the correct subplot.

a few things I added :

• apply a window to data before doing the fft (as the buffer of data does not start neither ends with zero)
• optional : use a exponential averaging of the fft data (a factor)

hope it helps

clc
load('ch1_first_5s.mat');
load('ch2_first_5s.mat');
fs = 128;
n1 = 0:1:((5 * fs) - 1);
n2 = 1:1:(fs - 1);

%%% FFT exponential averaging coefficient
a = 0 ;  % a = 0 : no averaging; a = 0.9 : highly averaged ; do not exceed 0.99 !!

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % my test data % debug only 
% dt = 1/fs;
% time = n1*dt;
% ch1_first_5s = 10*sin(2*pi*10*time);
% ch2_first_5s = 20*sin(2*pi*20*time);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


% data initialisation % 
ci = 1;
x_buffer = (ci + n2);
ch1_buffer = ch1_first_5s(ci + n2);
ch2_buffer = ch2_first_5s(ci + n2);
ch_buffer = [ch1_buffer;ch2_buffer];


%% plot 
h = figure(1),
subplot(2, 2, 1);
plot(x_buffer, ch_buffer(1,:), 'b');
title('Channel 1 - Time Domain EEG plot');
ylim([-200 200]); ylabel('Amplitude');
xlim([ci (ci + 128)]); xlabel('Discrete Time n (Sample Number)');
grid on;
subplot(2, 2, 3);
plot(x_buffer, ch_buffer(2,:), 'b');
title('Channel 2 - Time Domain EEG plot');
ylim([-200 200]); ylabel('Amplitude');
xlim([ci (ci + 128)]); xlabel('Discrete Time n (Sample Number)');
grid on;

N = length(n2)+1;
N1 = (2)^nextpow2(N);
window = hanning(length(n2));
X_k1 = fft(ch1_buffer.*window', N1); X_k1 = X_k1(1:(N1 / 2));
X_k2 = fft(ch2_buffer.*window', N1); X_k2 = X_k2(1:(N1 / 2));
Mag1_dB = 20*log10(abs(X_k1)*4/N1);
Mag2_dB = 20*log10(abs(X_k2)*4/N1);
Mag1_dB_old = Mag1_dB;
Mag2_dB_old = Mag2_dB;
freq = fs * (0:(N1 / 2) - 1) / N1; 
subplot(2, 2, 2);
plot(freq, Mag1_dB, 'b');
ylim([-60 40]); ylabel('Magnitude (dB)');
xlim([0 freq(end)]); xlabel('Frequency (Hz)');
title('Channel 1 - Freq Domain EEG plot');
grid on;
subplot(2, 2, 4);
plot(freq, Mag2_dB, 'b'); 
ylim([-60 40]); ylabel('Magnitude (dB)');
xlim([0 freq(end)]); xlabel('Frequency (Hz)');
title('Channel 2 - Freq Domain EEG plot');
grid on;

%%------The axis handle code-----------------------
h=findobj(gcf,'type','axes');

% h =   4×1 Axes array:
%   Axes    (Channel 2 - Freq Domain EEG plot)
%   Axes    (Channel 1 - Freq Domain EEG plot)
%   Axes    (Channel 2 - Time Domain EEG plot)
%   Axes    (Channel 1 - Time Domain EEG plot)
% please not that the axes array do not reflect the subplot order ! 
  
%% loop to update figure / axes handle
for ci = 1:1:511
    
    % update for LHS time plots
    x_buffer = (ci + n2);
    ch1_buffer = ch1_first_5s(ci + n2);
    ch2_buffer = ch2_first_5s(ci + n2);
    ch_buffer = [ch2_buffer;ch1_buffer]; % reversed order because so are the axes handle

   % update for RHS FFT plots
    X_k1 = fft(ch1_buffer.*window', N1); X_k1 = X_k1(1:(N1 / 2));
    X_k2 = fft(ch2_buffer.*window', N1); X_k2 = X_k2(1:(N1 / 2));
    Mag1_dB = a*Mag1_dB_old+(1-a)*20*log10(abs(X_k1)*4/N1); % exponential averaging 
    Mag2_dB = a*Mag2_dB_old+(1-a)*20*log10(abs(X_k2)*4/N1); % exponential averaging
    Mag_buffer = [Mag2_dB;Mag1_dB]; % reversed order because so are the axes handle
    % update FFT Mag 
    Mag1_dB_old = Mag1_dB;
    Mag2_dB_old = Mag2_dB;
    
%   update handles
    for k=1:4
      f=get(h(k),'children');
      if k == 3 || k == 4
          set(f,'xdata',x_buffer);
          set(h(k),'XLim',[ci (ci + 128)]);
          set(f,'ydata',ch_buffer(k-2,:));
      elseif k == 1 || k == 2
          set(f,'ydata',Mag_buffer(k,:));
      end
    end
     pause(0.05);
    drawnow
    
end

[Joshua Scicluna]

Hi, Thanks a lot for the above!!, I will take a look and let you know. Thanks