Extract signal components from mixed signal

30 次查看（过去 30 天）
Sohel Rana
Sohel Rana 2020-12-14
回答： Wayne King 2020-12-23
Hi,
I have signal that is consists of three individual signals. How can I separate these individual signals from the mixed signal. I did the FFT and then IFFT to get back to the original signal. The same way I wanna get the three individual signals (I_first, I_second, I_third) by performing FFT and IFFT.
The code I used:
clc
clear all;
close all;
format long
m=1000; I1=0.5; I2=0.3; I3=0.2; L1=100*m; L2=1000*m; n1=1; n2=1.446;
lam1=1530; lam2=1565;
inc=(lam2-lam1)/(2^10-1);
lam=lam1:inc:lam2;
Q12=(4*pi*n1*L1)./lam;
Q23=(4*pi*n2*L2)./lam;
Q13=Q12+Q23;
I_first=I1+I2+2*sqrt(I1*I2).*cos(Q12); % first signal
I_second=I2+I3+2*sqrt(I2*I3).*cos(Q23); % second signal
I_third=I1+I3+2*sqrt(I1*I3).*cos(Q13); % third signal
I=I1+I2+I3+2*sqrt(I1*I2).*cos(Q12)+2*sqrt(I2*I3).*cos(Q23)+2*sqrt(I1*I3).*cos(Q13); % Mixed signal
N=length(lam);
fs=1/inc;
dt=1/fs;
df=1/(N*dt);
f=(-N/2:N/2-1)*df;
xxx=5;
subplot(xxx,1,1)
plot(lam,I)
title('Mixed signal')
y=fft(I);
y1=fftshift(y);
y2=abs(y1);
subplot(xxx,1,2)
plot(f,abs(y))
title('FFT')
subplot(xxx,1,3)
plot(f,y2)
title('FFTshift')
y3=ifft(ifftshift(y1));
subplot(xxx,1,4)
plot(lam,abs(y3))
title('IFFTshift')

请先登录，再回答此问题。

采纳的回答

Wayne King
Wayne King 2020-12-23
编辑：Wayne King 2020-12-23
Hi Sohel, May I first suggest that you clean this code up a bit. It isn't easy to follow exactly what kind of signals you are creating. For example, instead of
lam1=1530; lam2=1565;
inc=(lam2-lam1)/(2^10-1);
Why not just ?
t = linspace(1530,1565,1024);
It appears from your code above you are using "lam" as your time vector. Then, you seemingly create a sine wave wtih expressions like in the case of I_first where I have used t for your "lam". What exactly are you trying to do here? A little explanation about what signal you are trying to create would help us determine if perhaps you have an inadvertent error in the signal model.
Note the equivalence of the following:
t = linspace(1530,1565,1024);
angle = 4*pi*1e5./t;
plot([angle' Q12'])
% Then you basically do
plot(cos(angle))
Were you trying to create a frequency-modulated signal with the above?
At any rate, I don't see how you can cleanly separate I_second from I_third. The frequencies of these two components are so close, whether intentional or not, that you have essentially just created one amplitude modulated signal. Note
plot([abs(fft(I_second))' abs(fft(I_third))'])
Now if I plot the sum of those two in time, you see the amplitude modulated signal
plot(I_second+I_third)
With respect to bandpass filtering, you are not going to be able to separate these components in a way that when you sum them back you get the original signal. You can however accomplish that with multiresolution techniques. In this case I would recommend a wavelet packet technique, modwptdetails. Please see Practical MRA for a introduction.
So here:
mra = modwptdetails(I,5,'fk18');
% first component
plot(mra(1,:))
% second and third component together
plot(mra(3,:))
Again, you cannot expect to separate I_second and I_third with any technique I know of. Now note that if I sum all the mra components back together
ts = sum(mra);
max(abs(ts-I))
I get back the original signal perfectly. Now if you compare the extracted MRA components, you see that except for the expected DC shifts (shifts in the mean), they quite accurately reproduce, I_first and I_second+I_third
subplot(2,1,1)
plot([mra(1,:)' I_first']), title('First');
axis tight
subplot(2,1,2)
plot([mra(3,:)' (I_second+I_third)']), title('Second+Third')
axis tight
To see that more clearly, let's add the DC shift in and compare the AM component extracted by the wavelet packet MRA with the original.
mu = mean(I_second+I_third);
figure
plot([mra(3,:)'+mu (I_second+I_third)']);
title('Comparison of wavelet packet with original')
axis tight
Hope that helps,
Wayne
  1 个评论
Sohel Rana
Sohel Rana 2020-12-23
Hi Wayne,
Thank you for your nice explanation and suggestion. My main goal is to separate I_first and I_second signals from the mixed signal (I) and then reconstrcut the orginal I_first and I_second. Yes, you are right, the frequencies of I_second and I_third are very close to each other but not same. When I look at the FFT of I_second and I_third, I can see their peaks are at different positions although close. I have only the parameters L1 and L2 that I can change willingly. A lot of works have been done on similar type of works. They used different types of digial filter to reconstrcut the original signal although their frequencies were very close to each other just like mine one. I beleive there should be some way to separate the I_first and I_second signals. I will look forward to getting your response.
https://ieeexplore.ieee.org/document/7997885
https://www.spiedigitallibrary.org/journals/optical-engineering/volume-53/issue-06/067101/Simultaneous-measurement-of-temperature-and-pressure-with-cascaded-extrinsic-FabryPerot/10.1117/1.OE.53.6.067101.full?SSO=1
https://ieeexplore.ieee.org/document/8558685

请先登录，再进行评论。

更多回答（2 个）

Abhishek Gupta
Abhishek Gupta 2020-12-17
Hi,
As per my understanding, you want to extract different components of a mixed signal. This task can be done in the following steps: -
  1. Retrieve frequencies of the components using Fast Fourier Transform (FFT)
  2. Use bandpass() or highpass() filter to extract different components out of the mixed signal.
Also, referring to the following MATLAB Answers, which might help you in achieving the task: -
  1 个评论
Sohel Rana
Sohel Rana 2020-12-17
编辑：Sohel Rana 2020-12-17
Hi Abhishek,
Thank you for your comments and suggestion. I have tried both ways but did not get the original signal. The reconstructed signals are quite different from the original one. Please see figure 1 (original) and figure 3 (reconstructed) in the code.
Code:
clc
clear all;
close all;
format long
m=1000; I1=0.5; I2=0.3; I3=0.2; L1=100*m; L2=1000*m; n1=1; n2=1.446;
lam1=1530; lam2=1565;
inc=(lam2-lam1)/(2^10-1);
%inc=0.1;
lam=lam1:inc:lam2;
Q12=(4*pi*n1*L1)./lam;
Q23=(4*pi*n2*L2)./lam;
Q13=Q12+Q23;
I_first=I1+I2+2*sqrt(I1*I2).*cos(Q12); % first signal
I_second=I2+I3+2*sqrt(I2*I3).*cos(Q23); % second signal
I_third=I1+I3+2*sqrt(I1*I3).*cos(Q13); % third signal
I=I1+I2+I3+2*sqrt(I1*I2).*cos(Q12)+2*sqrt(I2*I3).*cos(Q23)+2*sqrt(I1*I3).*cos(Q13); % Mixed signal
figure(1)
subplot(4,1,1);plot(lam,I_first); subplot(4,1,2);plot(lam,I_second);subplot(4,1,3);plot(lam,I_third);subplot(4,1,4);plot(lam,I);
%FFT
fs=1000;fn=fs/2;N=length(I);nfft=N;
FT=fft(I,nfft);
f=fs*(0:nfft-1)/nfft;
ff=fs*(0:nfft/2-1)/nfft;
FT2=FT(1:nfft/2);
[p,l]=findpeaks(abs(FT2),ff);
figure(2)
subplot(2,1,1);plot(f,abs(FT)); subplot(2,1,2);plot(ff,abs(FT2));
% filtering
ps=2; kk=l;k=kk(ps);dn=0.0002; f1=k-dn; f2=k+dn;
order=50; cutoff=0.01;
z3=highpass (I,42,fs); % since l(3)= 43.945312500000000
z2=bandpass(I,[f1 f2],fs);
z1=I-(z2+z3);
Z=z1+z2+z3;
figure(3)
subplot(4,1,1);plot(lam,z1); subplot(4,1,2);plot(lam,z2);subplot(4,1,3);plot(lam,z3);subplot(4,1,4);plot(lam,Z);

请先登录，再进行评论。

Wayne King
Wayne King 2020-12-23
Hi Sohel, for the given frequencies and signal lengths, the spectra of I_second and I_third are not delta functions in frequency. So the fact that their peaks are slightly different isn't going to be sufficient. If you knew they were sine waves and knew what their frequencies were, perhaps it would be possible to resample them in such as way that their respective frequencies (with all their associated energy) fell exactly on a DFT bin (and separate DFT bins) and then you could perhaps separate them, but here there is considerable overlap so I don't see how you are going to separate them as you hope to. I'm sorry I don't have time to read those papers, but if those papers contain the work of people constructing signals exactly like yours and claiming to separate I_second and I_third perfectly, perhaps contacting them directly? Maybe they would help.
Wayne

类别

在 Help Center 和 File Exchange 中查找有关 Image Category Classification 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by