- 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)
Animated magnitude spectrum (windowed fft)
3 次查看(过去 30 天)
显示 更早的评论
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;
0 个评论
采纳的回答
Mathieu NOE
2021-2-8
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 :
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
更多回答(0 个)
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 EEG/MEG/ECoG 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!