Injecting (and combining) a sine wave after a period of time

Hello everyone,

So,I have this one signal of 100s seconds long. It is an EEG signal which looks very much like a random signal. I want to inject a sine wave tACS data (5s long sine wave) into it at 30s. But when I say injection, I didn't mean putting sine signal over another EEG signal. At the time 30-35 seconds, I wanted to combine the two signals via addition. Lastly, displaying the full 100s EEG signal with 30-35s: sine wave data (tACS)+EEG. I have attached few image to pictorially what I meant. Here is the current code, that I am working on:

load('Test_Sham_Alpha.mat');%load EEG data file
%plot(EEG.Time(1,1:46600)' , EEG.Data(1,1:46600)')

%define sine wave
%%Time specifications:
   Fs = 500;                   % samples per second
   dt = 1/Fs;                   % seconds per sample
   StopTime = 5;             % seconds
   t = (0:dt:StopTime-dt)';     % seconds
   %%Sine wave:
   Fc = 5;                     % hertz
   A= 0.001;
   x = A*sin(2*pi*Fc*t);
   % Plot the signal versus time:
   %figure;
   %plot(t,x);%length t elements = 2500
   %xlabel('time (in seconds)');
   %title('tACS Signal');

%inject f 5 Hz tACS at time: 30 seconds of EEG data (Y= 30453934) 
%for 5 seconds (25 cycles)
time= EEG.Time(1,1:46600)
for EEG.Time(1,15000:17500)%length of EEG.time =2501
    %add y-axis value for both signals 
   y1= EEG.Data(1,1:46600);
   y2= x
   y= y1+y2
end

%i am not sure if the above make sense or now
%not sure how to continue 

%for the figure: I wanted that the plot shows EEG of its full length data
%where between 30-35s of EEG, there is a signal combination (adding of tACS sine
%wave with the EEG)and the rest of EEG signal remained as it is
figure ('Name','Combined EEG+ tACS');
plot(t,y);%length t elements = 2500
xlabel('time (in seconds)');

I have also attached the EEG data that I extracted into an excel file (but only from 25s -40s because 100s is too big but it can easily be plot onto Matlab (fast too), but it is data that requires permission before posting it on public), if you have an idea on how to do it. There are other things that I want to add it on data: ramp-up features, and implementing envelope onto the signals etc. But before I go on to those things, I'd like to settle the basic issues. Note: I can do this on Excel but it couldn't handle that too much data, so I wanted to use maltab but I am not very good at the coding part. I am pretty sure it should be simple. I thought of cutting the datas and then concatenate them but I don't think that is efficient.

Please asisst me. I reallyyyy need help on this. I feel losst :'(

Regards,

Anis

Answer 1

[Rohit Pappu]

A plausible solution is as follows

• Generate the required sinusoidal signal

%define sine wave
%%Time specifications:
   Fs = 500;                   % samples per second
   dt = 1/Fs;                   % seconds per sample
   StopTime = 5;             % seconds
   t = (0:dt:StopTime-dt);     % seconds , assuming t is a 1xn vector to match the dimensions of EEG data
   %%Sine wave:
   Fc = 5;                     % hertz
   A= 0.001;
   x = A*sin(2*pi*Fc*t);
   

• Select the EEG data corresponding to the time stamps required (In this case, index 15001 to 17500)

y_to_be_modified = EEG.Data(1,15001:17500);

• Add the sine wave to the EEG data

y_injected = x+y_to_be_modified;

• Replace the original EEG data section with this injected data

EEG.Data(1,15001:17500) = y_injected;

• Plot the EEG Data

plot(EEG.Time(1,1:46600) , EEG.Data(1,1:46600))

[Anis Fatini Abdul Aziz]

I tried this but it does not work. The sine wave was not added onto the 5s window so when I plot them they look the same (no difference). Here is the code I used:

%define sine wave
%%Time specifications:
 Fs = 500;                   % samples per second
 dt = 1/Fs;                   % seconds per sample
 StopTime = 5;             % seconds
 t = (0:dt:StopTime-dt);     % seconds , assuming t is a 1xn vector to match the dimensions of EEG data
 Fc = 5;                     % hertz
 A= 0.001;
 x = A*sin(2*pi*Fc*t);
   
load('Test_Sham_Alpha.mat');%load EEG data file

y_to_be_modified = EEG.Data(1,15001:17500);
y_new= y_to_be_modified*10^-9
y_injected = x+y_new;
EEG.Data(1,15001:17500) = y_new;
figure(2);
plot(EEG.Time(1,15001:17500) , y_new)

I have attached the data file with it: Note that the y-axis should be in nV (E10-9) so hence that part of the code. I am suspecting it because of the sine wave is not in vector form? Like in data set? But I am not sure, if so, any suggestion? Should I make it into a MATLAB data first like put them into tables so the addition makes more sense?

[Rohit Pappu]

There are a few issues with the above code

• Line 18 :Only a sample of the signal is being plotted, not the entire signal. To plot the entire signal

plot(EEG.Time(1:46600), EEG.Data(1:46600));

• Line 14: Only the injected part is converted to nV, not the entire EEG Data. This can lead to skewed graph like this

An alternative would be as follows

%% Scale the entire EEG.Data to the same units as x
EEG.Data = scalingFactor*EEG.Data;
y_to_be_modified = EEG.Data(1,15001:17500);
y_injected = x+y_to_be_modified; % Assuming x and y_to_be_modified have the same units ie. nV
EEG.Data(1,15001:17500) = y_injected;
figure(3);
plot(EEG.Time(1,1:46600) , EEG.Data(1,1:46600)) % Plots the entire signal

[Anis Fatini Abdul Aziz]

Hi,

Thank you very much for your clarification and assistance! I really appreciate it. It makes a lot more sense. It works (attached file: test_matsin.m, same code as above). The plot looks similar as a reference file I used (attached file: eeg+tacs_ref). However, this is a case whereby I am taking a matlab-generated sine wave. I was only doing this for comparison and testing purposes.

In my case, what I actually want to do is to use a sine wave result file that I got from a simulation (attached file: mymat.mat). I got to do that already, the code only differs slightly (where the sine wave that was previously generated, replaced with the upper part shown):

load('mymat.mat');
mili =10^-3
tacs= mili*(q(1:30000,1)');

load('eeg_sham.mat');%load EEG data file
scalingFactor = 10^-9
%% Scale the entire EEG.Data to the same units as x
EEG.Data = scalingFactor*EEG.Data;
y_to_be_modified = EEG.Data(1,15001:45000);
y_injected = tacs+y_to_be_modified; % Assuming x and y_to_be_modified have the same units ie. nV
EEG.Data(1,15001:45000) = y_injected;
figure(3);
plot(EEG.Time(1,1:46600) , EEG.Data(1,1:46600)) % Plots the entire signal

So, right now, I want to extract the envelope of the reference signal and apply it on to the signal that I have plotted, i.e. the output signal of code above. Is that possible? I did the envelope part already and I think output is partly odd (where is the transients), I am expecting something like Figure 1(attached). Now, I do not know how to apply the extracted envelope onto the signal. Any ideas?

load('eeg+tacs_ref_5Hz_1mA.mat');%load EEG data file
figure(1);
plot(EEG.Time(1,1:50400) , EEG.Data(1,1:50400))
t= EEG.Time(1,15001:16201) % test smaller window size (did a big one, not possible)
x= EEG.Data(1,15001:16201)

%for envelope extraction test 
y = hilbert(x);
env = abs(y);
plot_param = {'Color', [0.6 0.1 0.2],'Linewidth',2}; 

figure (2)
plot(t,x)
hold on
plot(t,[-1;1]*env,plot_param{:})
hold off
title('Hilbert Envelope')

%y-axis is readily scaled to Volts
load('mymat.mat');
mili =10^-3
tacs= mili*(q(1:30000,1)');
   
load('Test_Sham_Alpha.mat');%load EEG data file
scalingFactor = 10^-9
%% Scale the entire EEG.Data to the same units as x
EEG.Data = scalingFactor*EEG.Data;
y_to_be_modified = EEG.Data(1,15001:45000);
y_injected = tacs+y_to_be_modified; % Assuming x and y_to_be_modified have the same units ie. nV
EEG.Data(1,15001:45000) = y_injected;
figure(3);
plot(EEG.Time(1,1:46600) , EEG.Data(1,1:46600)) % Plots the entire signal, let's call it channel 1

%appply envelope to signal on channel 1
% add? multiply? insert window trace? 

Please assist me. Thank you very much in advance!

[Rohit Pappu]

For obtaining the envelope of the target signal, Signal Envelope can be used.

In above case:

[y_upper, y_lower] = envelope(y_injected) % Generate the envelope for injected values
plot(EEG.Time(1,1:46600) , EEG.Data(1,1:46600)) % Plot the entire signal
hold on
plot(EEG.Time(1,15001:45000),  y_upper, y_lower) % Overlay the envelope on the previous plot
hold off

[Anis Fatini Abdul Aziz]

Hi Rohit,

Thank you so much for the prompt reply. I tried that but I got an error shown below, whereby the plot function cannot be used with two y-values?

Also, from my last post, I attached an image of expected envelope: Figure 1 file. And then I ask, why does it look like that? As in why does it not detect the transient or the ramp up part of the envelope. When I use a bigger window where I include the before injection (before 30 s) till 40s, I can see the ramp up part now (Figure 2 attached).

I have got a que: How do I plot the envelope on original signal, in the same figure? I tried doing that (code below), but it came out really odd (Figure 3 attached)

load('Test_1mA_5Hz_Alpha.mat');%load EEG data file
figure(1);
plot(EEG.Time(1,1:50400) , EEG.Data(1,1:50400))
t= EEG.Time(1,14001:20001) %use a bigger window size 
x= EEG.Data(1,14001:20001)

%for envelope extraction test 
y = hilbert(x);
env = abs(y);
plot_param = {'Color', [0.6 0.1 0.2],'Linewidth',2}; 

figure (2)
plot(EEG.Time(1,1:50400) , EEG.Data(1,1:50400))
hold on 
plot(t,x)
hold on
plot(t,[-1;1]*env,plot_param{:})
hold off
title('Hilbert Envelope')

%y-axis is readily scaled to Volts
load('mymat.mat');
mili =10^-3
tacs= mili*(q(1:30000,1)');
   
load('Test_Sham_Alpha.mat');%load EEG data file
scalingFactor = 10^-9
%% Scale the entire EEG.Data to the same units as x
EEG.Data = scalingFactor*EEG.Data;
y_to_be_modified = EEG.Data(1,15001:45000);
y_injected = tacs+y_to_be_modified; % Assuming x and y_to_be_modified have the same units ie. nV
EEG.Data(1,15001:45000) = y_injected;
figure(3);
plot(EEG.Time(1,1:46600) , EEG.Data(1,1:46600)) % Plots the entire signal, let's call it channel 1

%appply envelope to signal on channel 1
figure(4)
[y_upper, y_lower] = envelope(y_injected) % Generate the envelope for injected values
plot(EEG.Time(1,1:46600) , EEG.Data(1,1:46600)) % Plot the entire signal
hold on
plot(EEG.Time(1,15001:45000),  envelope(y_injected)) % Overlay the envelope on the previous plot
hold off

Thank you so much. I really appreciate your help.

Best Regards,

Anis

[Rohit Pappu]

For plotting multiple signals on the same set of axes, please use hold. To customise the line colors , refer to this documentation.

A brief outline would be as follows:

[env_upper,env_lower] = envelope(EEG.Data(1,1:46600))
plot(EEG.Time(1,1:46600) , EEG.Data(1,1:46600)) % Plot the entire signal
hold on
plot(EEG.Time(1,1:46600), env_upper, EEG.Time(1,1:46600), env_lower) % Overlay the envelope on the above plot. Use line spec to change the line colors
hold off

[y_upper, y_lower] = envelope(y_injected) % Generate the envelope for injected values
hold on
plot(EEG.Time(1,15001:45000),  y_upper, EEG.Time(1,15001:45000), y_lower) % Overlay the envelope of injected
% signal on the previous plot . Use Line spec to change the color of the line
hold off

Additional documentaion can be found here

[Anis Fatini Abdul Aziz]

Hi Rohit,

This only overlays the envelope ontop of the signal right? What I am trying to do is to apply the envelope onto the signal i.e. get the signal to follow the same envelope shape/ pattern as the reference signal?

Regards,

Anis

[Anis Fatini Abdul Aziz]

I wrote the code again but hmm, it is only overlaying and it looks messed up (suddenly it is at the bottom):

%y-axis is readily scaled to Volts
load('channel_1_F4.mat');
mili =10^-3
tacs= mili*(q(1:30000,1)');
   
load('Test_Sham_Alpha.mat');%load EEG data file
scalingFactor = 10^-9
%% Scale the entire EEG.Data to the same units as x
t_signal_A =EEG.Time(1,15001:45000)
EEG.Data = scalingFactor*EEG.Data;
y_to_be_modified = EEG.Data(1,15001:45000);
y_injected = tacs+y_to_be_modified; % Assuming x and y_to_be_modified have the same units ie. nV
EEG.Data(1,15001:45000) = y_injected;
figure(1);
plot(EEG.Time(1,1:46600) , EEG.Data(1,1:46600)) % Plots the entire signal, let's call it channel 1

load('eeg+tacs_ref_5Hz_1mA.mat');%load reference signal for envelope
ref_t=EEG.Time(1,1:50400)
ref_y = EEG.Data(1,1:50400)
figure(2)
[env_upper,env_lower] = envelope(ref_y)
plot(ref_t , ref_y) % Plot the entire signal
hold on
plot(ref_t, env_upper, ref_t, env_lower) % Overlay the envelope on the above plot. Use line spec to change the line colors
hold off
[y_upper, y_lower] = envelope(y_injected) % Generate the envelope for injected values
hold on
plot(t_signal_A,  y_upper, t_signal_A, y_lower) % Overlay the envelope of injected
% signal on the previous plot . Use Line spec to change the color of the line
hold off

[Rohit Pappu]

Please refer to this documentation - Envelope extraction

[Anis Fatini Abdul Aziz]

I have gone through the documentation. Is it not the same as what I have done before? I don't see how this helps to make a certain signal follow another's signal's envelope. Please clarify this matter.

[Rohit Pappu]

Could you please explain, what does making a certain signal follow another signal's envelope mean?

Thank you

[Anis Fatini Abdul Aziz]

I have two signal: signal A (modified signal) and signal B (eeg+tacs_ref_5Hz_1mA.mat). I extract the envelope of signal B. My aim is to have signal A takes (follows/ trace) the pattern/ shape of signal B's extracted envelope. Also, actually, for the extracting envelope of signal B (call it X), the important part to extract/trace is between 30s to 90s, the rest of signal B is not needed. In the end, I am expecting the window time from 30s-90s of signal A follows/have same pattern of X. Does that make sense?

[Rohit Pappu]

Based on the above problem statement, a possible solution is as follows

• B_extract = B(30:90)
• Modulate Signal B with A ie. B_extract+A.*B_extract (Amplitude modulation)

[Anis Fatini Abdul Aziz]

Thank you for your idea. I see, alright I will try that and update again. I just have one question, that would that signal A, B and B_extract has to be in the same units right? If not, I think when I modulate its amplitude and plot the entire thing it would be weird, isn't it?

[Rohit Pappu]

Yes, they all need to be in the same unit

[Anis Fatini Abdul Aziz]

Wait, I don't understand something. For this one:

• Modulate Signal B with A ie. B_extract+A.*B_extract (Amplitude modulation)

The B extract is the envelope right? The envelope has upper and lower part. So how do I get the signal A's upper and lower part separately before I can do the amplitude modulation using the formula you gave. Or did I understand this the wrong way? Because I don't think I can straight away mutiply the extract with signal? because the envelope has the two part?

Also, shouldn't it be modulate signal A with envelope B? beccause that is what I am trying to do? modify signal A to follow same pattern as B. So B is like the carrier signal? Right?

[Anis Fatini Abdul Aziz]

Okay so I tried something, but it looks so weird, like there is AM on the top part of the signal. Also, why are there no ramp up as I expect there would be? Idk T.T Please help, where did I do wrong. Here is my code:

load('Test_1mA_5Hz_Alpha.mat');%load EEG data file
scalingFactor = 10^-9
env_window = scalingFactor* EEG.Data(1,15001:45000)%window of interest for envelope extraction 
t_window = EEG.Time(1,15001:45000)% 30 s to 90 s
t_ref = EEG.Time(1,1:50400)%entire signal from 0s to 100s
y_ref = scalingFactor*EEG.Data(1,1:50400) 
figure (1);
plot( t_ref,y_ref ) 
title('Reference: EEG + 1mA  5Hz tACS Signal from Phantom Head')
xlabel ('Time, s');
ylabel ('Voltage, V');

%Using ROAST tACS voltage results recorded
%y-values of tacs is scaled from mV to V and transposed
load('channel_3_F3.mat');
mili =10^-3
tacs= mili*(q(1:30000,1)');
   
y_injected = tacs+y_to_be_modified; % Assuming x and y_to_be_modified have the same units ie. nV
EEG.Data(1,15001:45000) = y_injected;
figure(2);
plot(EEG.Time(1,1:46600) , EEG.Data(1,1:46600)) % Plots the entire signal
title('Combined EEG+ROAST-tACS')
xlabel ('Time,s');
ylabel ('Voltage, V');

figure(3)
[env_upper,env_lower] = envelope(env_window)
plot(t_window , env_window) % Plot the entire signal
hold on
plot(t_window, env_upper, t_window, env_lower) % Overlay the envelope on the above plot. Use line spec to change the line colors
hold off
title('Envelope of Referred Signal from Phantom Head')
xlabel ('Time,s');
ylabel ('Voltage, V');

%possible idea
%B_extract = B(30:90)
%Modulate Signal B with A ie. B_extract+A.*B_extract (Amplitude modulation)
load('Test_Sham_Alpha.mat');
EEG.Data = scalingFactor*EEG.Data;
am_mod = envelope(env_window) + y_to_be_modified.*envelope(env_window);
EEG.Data(1,15001:45000) = am_mod;
figure (4)
plot(EEG.Time(1,1:46600) , EEG.Data(1,1:46600))
title('Amplitude Modulated EEG+tACS Signal')
xlabel ('Time,s');
ylabel ('Voltage, V');

Please assist :'(

[Anis Fatini Abdul Aziz]

does anyone know how?? please, it is still not working.