Subtracting complex number from a discrete signal

Hello Matlab community,

I'm trying to remove a background signal from a set of discrete signals. I have determined that the background signal is represented by -0.000053 + 0.0030995i, or 0.0031*e^(i*1.588). In my code, the "j" loop represents each of the 11 positions where signals were recorded. I'm having trouble understanding how to actually perform the subtraction of the background signal from the 11 discrete signals; it seems that I need to convert the discrete digital signals into the complex plane but I don't know how to do that. Any assistance would be appreciated; I've attached the code.

Thank you,

Michelle

% DataF.mat has 250000 data points per position.  Each point of data
% corresponds to 1mm of motion beginning at -5mm from center to 5mm from 
% center ranging 10mm total, resulting in 11 data points.
%
clc;
clear;
close all;

%% Constants
EPSfreq=2.36e3;             %EPS carrier                      
rate=50000;                 %LabView data file Sample Clock rate = rate 
dt=1/rate;
points=250000;              

%% Main Code
load('DataF.mat');
% Initialized Matrices
% [c,t,tm,out]=deal(zeros(points,numPoints));

% nSamp=Temp_points
for i=1:points
    t(i)=i*dt;
    tm(i)=t(i)*1e6;                                                        %Time in microseconds
end
%Position loop: 11 positions, 100 points per measured location
for j=1:11                                                                 
%Output from filter (pickup signal)
out(:,j) = DataF((j-1)*points+1:(j*(points)), 1);                      
%Reference from excitation coils
c(:,j) = DataF((j-1)*points+1:(j*(points)), 2);        

%Function detrend removes data bias by computing mean values and subtrancting them from signals
%Use if DAQ system puts offset on signals
out_d(:,j) = detrend(out(:,j));                                                       
c_d(:,j) = detrend(c(:,j));   

%Plot of detrended signals: plots for each position
% figure
% plot(t,c_d(:,j),t,out_d(:,j)); 
% xlabel('Time (s)')
% ylabel('Amplitude (V)')

%Cross-correlation and auto-correlation
    %Auto correlation reference
    Rx(j) = (1/points)*sum(c_d(:,j).^2);                                                   
    %Auto correlation measured                                                      
    Ry(j) = (1/points)*sum(out_d(:,j).^2);     

    %Amplitude of reference
    Ax(j) = sqrt(2*Rx(j)); 
    %Amplitude of measured
    Ay(j) = sqrt(2*Ry(j));                                                 
    
    %Cross-correlation             
    R(j) = (1/points)*sum(c_d(:,j).*out_d(:,j));    
    %Relative phase in degrees
    relphase(j)=acos(R(j)/sqrt(Rx(j)*Ry(j)))*(180/pi);     

    %Amplitude ratio
    relamp(j) = Ay(j)/Ax(j);                                               
    end
%Sign adjust loop
for k=1:5                                                                 %First 5 positions
    signadj(k)=sign(relphase(k)-180);                                      %Adjusts amplitude graph slope
    relamp(k)=relamp(k)*signadj(k);                                        %Output results in linear graph slope
end

figure
plot([-5:1:5],relamp,'-*');grid
xlabel('Displacement')
ylabel('Derived relative amplitude')
% axis([-5 5 -0.05 0.05]) 
figure
plot([-5:1:5],relphase,'-*');grid
xlabel('Displacement')
ylabel('Phase difference')
figure
subplot(2,2,1);
plot(Rx)
title('Auto Correlation Sum: Excitation')
legend('X')
grid on
subplot(2,2,2);
hold on
plot(R)
plot(Ry)
title('Auto and Cross Correlation Sum: Cross correlation and Pickup')
legend('R','Y')
figure
subplot(2,2,1);
plot(Ax)
title('Excitation Amplitude')
legend('Ax')
grid on
subplot(2,2,2);
plot(Ay)
title('Pickup Amplitude')
legend('Ay')
grid on

[Mathieu NOE]

hello Michelle

can you share the data (mat file) as well ?

tx

[Michelle Weinmann]

My data file is 30 MB, even after compressing in a zip file. Matlab won't let me attach it since it's larger than 5 MB. Each signal has 250,000 samples for each 11 data points. Is there a way to attach large files?

[gotjen]

Can you crop out just one of the discrete signals and upload that?

[Michelle Weinmann]

Good idea, I've attached the first discrete signal. The first column is the signal named out and the second is the signal named c.

[Walter Roberson]

Are your signals real-valued or complex-valued? If they are real-valued, then the background noise is certain to be real-valued, not complex-valued. If they are complex-valued then just subtract the complex number from the array of signals.

[Michelle Weinmann]

My signals are real-valued. The background signal isn't noise on the signal, it's a separate signal with a different phase. The data I collected was from a sensor similar to an LVDT, I should have measured the null case but I didn't know to at the time; there is always an unrelated background signal it picks up from ambient electromagnetic sources that needs to be extracted for the analysis to be correct. Let me know if this wasn't helpful.

[Mathieu NOE]

the provided file contains only real valued data. No sign of complex background signal

find(abs(imag(data))> eps)
ans =
  0×1 empty double column vector

[Michelle Weinmann]

I might not have explained this correctly, forgive me I'm new to DSP. The background signal isn't complex, but I've been told by my advisor that in order to subtract one signal from another they both have to be represented in the complex plane. I've converted the background signal into a complex representation as Acos(wt+phi) = Ae^(j*theta), but the signal itself is real. I would like to subtract the background signal from my discrete data signal, but I don't know how to convert the discrete data signal into a complex representation and then perform the subtraction. The outputs of the script attached are the amplitude and phase at each location, so I can't use those to manipulate my discrete signals into complex representation. Please let me know if you have more questions. Thank you all for looking into my question.

[Mathieu NOE]

if you have electric perturbations from mains , this happen to known and fix frequencies (50 or 60 Hz and maybe a bit of harmonics)

why not simply use a notch filter and do the filtering directly on the time data ?

example attached

Answer 1

[Walter Roberson]

A = 123

A = 123

B =  -0.000053 + 0.0030995i

B = -0.0001 + 0.0031i

C = A - B

C = 1.2300e+02 - 3.0995e-03i

In other words, you do not need to do anything special to convert the signal into the complex plane: MATLAB will automatically extend real values to complex as needed.

[Walter Roberson]

You seem to be under the impression that you must do something specific to convert your discrete signal to complex representation. If so then

complex_signal = complex(real_signal, 0)

If you need complex form A*exp(j*theta) then this is the same as theta = 0, as clearly exp(j*0) --> exp(0) -> 1 and A*1 is A

However, as I pointed out, in MATLAB you can simply subtract the complex value and MATLAB will automatically extend the real value to complex .

Answer 2

[Michelle Weinmann]

Hello all,

It turns out I had the wrong idea about what I needed to do. I needed to subtract a complex number representing a signal from another complex number representing the output signal using the output values "relamp" as amplitude and "relphase" as phase. I have adjusted my code as required and attached it for reference.

However, it's not performing quite how it should. According to my advisor, the "relamp" values should be quite close to linear and I have redundant steps in the subtraction step which changes some of the values slightly; the linear relationship is slightly off. I can't see where the redundant steps are and how they would affect the relamp values. I've also attached a shortened data file (DataResized.mat) and the plot resulting (RelativeAmplitudePlot.fig); the first half of the points look ok but the second half do not. I've also attached a photo of what the plot should look like (NewPlot.jpg). Any help would be appreciated.

Thank you,

Michelle

% DataF.mat has 250000 data points per position.  Each point of data
% corresponds to 1mm of motion beginning at -5mm from center to 5mm from 
% center ranging 10mm total, resulting in 11 data points.
%
clc;
clear;
close all;

%% Constants
EPSfreq=2.36e3;             %EPS carrier                      
rate=50000;                 %LabView data file Sample Clock rate = rate 
dt=1/rate;
points=2500;              

%% Main Code
load('DataResized.mat');
% Initialized Matrices
% [c,t,tm,out]=deal(zeros(points,numPoints));

% nSamp=Temp_points
for i=1:points
    t(i)=i*dt;
    tm(i)=t(i)*1e6;                                                        %Time in microseconds
end
%Position loop: 11 positions, 100 points per measured location
for j=1:11                                                                 
%Output from filter (pickup signal)
out(:,j) = DataRsz((j-1)*points+1:(j*(points)), 1);                      
%Reference from excitation coils
c(:,j) = DataRsz((j-1)*points+1:(j*(points)), 2);        

%Function detrend removes data bias by computing mean values and subtrancting them from signals
%Use if DAQ system puts offset on signals
out_d(:,j) = detrend(out(:,j));                                                       
c_d(:,j) = detrend(c(:,j));   

%Plot of detrended signals: plots for each position
% figure
% plot(t,c_d(:,j),t,out_d(:,j)); 
% xlabel('Time (s)')
% ylabel('Amplitude (V)')

%Cross-correlation and auto-correlation
    %Auto correlation reference
    Rx(j) = (1/points)*sum(c_d(:,j).^2);                                                   
    %Auto correlation measured                                                      
    Ry(j) = (1/points)*sum(out_d(:,j).^2);     

    %Amplitude of reference
    Ax(j) = sqrt(2*Rx(j)); 
    %Amplitude of measured
    Ay(j) = sqrt(2*Ry(j));                                                 
    
    %Cross-correlation             
    R(j) = (1/points)*sum(c_d(:,j).*out_d(:,j));    
    %Correlation ratio for relphase
    Rrat(j) = R(j)/sqrt(Rx(j)*Ry(j));
    %Relative phase gives angles between 0 and 180 deg
    relphase(j)=acos(Rrat(j));     

    %Amplitude ratio
    relamp(j) = Ay(j)/Ax(j);  
    
    %Subtract background signal in complex plane
    V1(j)= relamp(j)*exp(1i*(relphase(j)));
    %Background signal
    V2 = -0.00005332875 + 0.0030995i;
    V3(j) = V1(j)-V2;
    %Update relamp and relphase values once background signal subtracted
    relamp(j) = abs(V3(j));
    relphase(j) = angle(V3(j))*180/pi;
    %Update pickup amplitude once background signal subtracted
    Ay(j) = relamp(j)*Ax(j);
    end
%Sign adjust loop: since acos only outputs angles from 0 to 180
for k=1:11                                                                 
    if Rrat(k) > 0
        relphase(k) = relphase(k) - 2*pi;
    end
    if Rrat(k) < 0
       relamp(k) = -relamp(k); 
    end
end

figure
plot([-5:1:5],relamp,'o-');grid
legend('Adjusted pickup relative amplitude','V3')
xlabel('Displacement')
ylabel('Derived relative amplitude')