How to apply Inverse FFT after filtering to get the original signal back?

92 views (last 30 days)
I have applied FFT on a mixed signal to find its component. Then I used low pass filtering to separate specifc signal component. However, when I applied inverse IFFT, I did not get the original signal back (the signal before the FFT). I think I made a mistake while applying inverse FFT. Could you please help me with that?
Here is 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=1520; lam2=1580;
inc=(lam2-lam1)/(2^12-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=500;fn=fs/2;N=length(I);nfft=N;
FT=fft(I);
f=fs*(0:nfft-1)/nfft;
ff=fs*(0:nfft/2-1)/nfft;
FT2=FT(1:nfft/2); % positive half
[p,l]=findpeaks(abs(FT2),ff);
figure(2)
subplot(2,1,1);plot(f,abs(FT)); subplot(2,1,2);plot(ff,abs(FT2));
%Hamming window
figure(3)
f1=l(1); f2=l(2); f3=l(3); freq=[f1 f2 f3];
mm=(0.1*f1)/(fs/2);%define tansition bandwidth
M=round(8/mm);%define the window length
NN=M-1;%define the order of filter
b=fir1(NN,1*f1/(fs/2));%use the firl function to design a filter
[h,f]=freqz(b,1,512);%amplitude-frequency characteristic graph
plot(f*fs/(2*pi),20*log10(abs(h)))
xlim([0 10]);
%wavelength domain
figure(4)
sf=filter(b,1,I);
subplot(3,1,1)
plot(lam,sf) % plot wavelength domain signal after filtering
Fsf=fft(sf,512);%frequency-domain diagram after filtering
AFsf=abs(Fsf);%the amplitude
f=(0:255)*fs/512;%frequency sampling
subplot(3,1,2)
plot(f,AFsf(1:256))%plot the frequency domain diagram after filtering
y=ifft(Fsf); % IFFT
subplot(3,1,3)
plot(y)

Accepted Answer

Star Strider
Star Strider on 19 Dec 2020
You would have to use fftshift and then apply the results of freqz symmetrically to each half of the complex fft result, so using it on the right half and flipping it for the left half, then doing element-wise multiplication, then using ifftshift to reverse it, and then doing the inverse fft. (I have done that in another Answer that I cannot now locate, and it is more trouble than it is worth.)
A much easier way to do what you want (with all the complicated bits already solved) is to use the fftfilt function. It will use your ‘b’ vector from the FIR filter you already designed. I strongly recommend that approach!
  14 Comments
Sohel Rana
Sohel Rana on 21 Dec 2020
Thank you. I'm not much concern about the amplitude. My main focus is to get the similar shape of the graph. The reconstructed graph shape should be similar to the I_second not the I.
While I'm putting the Wp = [1.216 1.217]/Fn; it's showing the following error.
Error using fft
FFT length must be a nonnegative integer scalar.
Error in freqz>firfreqz (line 210)
h = fft(b,s.*nfft).';
Error in freqz (line 127)
[h,w,options] = firfreqz(b/a,options);
Error in star22 (line 42)
freqz(sos, 2^20, Fs)
Star Strider
Star Strider on 21 Dec 2020
My code ran without error (in R2020b, Update 3), including the slight change in ‘Wp’. I have no idea what the problem could be. If there had been problems, I would have fixed them and noted any changes before I posted it.

Sign in to comment.

More Answers (2)

Siddharth. A
Siddharth. A on 1 Jan 2021
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=1520; lam2=1580; inc=(lam2-lam1)/(2^12-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=500;fn=fs/2;N=length(I);nfft=N; FT=fft(I); f=fs*(0:nfft-1)/nfft; ff=fs*(0:nfft/2-1)/nfft; FT2=FT(1:nfft/2); % positive half [p,l]=findpeaks(abs(FT2),ff); figure(2) subplot(2,1,1);plot(f,abs(FT)); subplot(2,1,2);plot(ff,abs(FT2));
%Hamming window figure(3) f1=l(1); f2=l(2); f3=l(3); freq=[f1 f2 f3]; mm=(0.1*f1)/(fs/2);%define tansition bandwidth M=round(8/mm);%define the window length NN=M-1;%define the order of filter b=fir1(NN,1*f1/(fs/2));%use the firl function to design a filter [h,f]=freqz(b,1,512);%amplitude-frequency characteristic graph plot(f*fs/(2*pi),20*log10(abs(h))) xlim([0 10]);
%wavelength domain figure(4) sf=filter(b,1,I); subplot(3,1,1) plot(lam,sf) % plot wavelength domain signal after filtering Fsf=fft(sf,512);%frequency-domain diagram after filtering AFsf=abs(Fsf);%the amplitude f=(0:255)*fs/512;%frequency sampling subplot(3,1,2) plot(f,AFsf(1:256))%plot the frequency domain diagram after filtering y=ifft(Fsf); % IFFT subplot(3,1,3) plot(y)

Sohel Rana
Sohel Rana on 4 Jan 2021
Hi Star,
My apologies for asking tons of questions.
I have applied the attached code to reconstruct the signal but could not get the exact signal. If you run the code, you will see the length before and after reconstuct are different (Fig. 5). If it is not same, at least it should be very close to the original length. Few things that I'm not understanding:
  1. How could we fix the passband and stopband ripple? The shape of the reconstructed signals change with these values.
  2. The range of lam can be varied at will. L1 and L2 can also be varied at will. Feel free to change.
  3. I'm not concern about the amplitude of the reconstcued signals. What I need is to get the similar shape so that the length of L1 and L2 become same both for original and reconstrcuted signals.
  4. Most of the paper I read did FFT and then used hamming window, digital filtering to get the signal back.
  5. Could you please tell me how can I get the reconstructed signals similar to the original signals?
I would be very grateful to you if you help me with this.
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=1520; lam2=1560;
inc=(lam2-lam1)/(2^15-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
%I=I1+I2+2*sqrt(I1*I2).*cos(Q12)+2*sqrt(I2*I3).*cos(Q23); % Mixed signal
figure(1)
subplot(4,1,1);plot(lam,I_first);title('First Signal'); subplot(4,1,2);plot(lam,I_second);title('Second Signal');subplot(4,1,3);plot(lam,I_third);title('Third Signal');subplot(4,1,4);plot(lam,I);title('Total Signal');
%subplot(3,1,1);plot(lam,I_first);title('First Signal'); subplot(3,1,2);plot(lam,I_second);title('Second Signal');subplot(3,1,3);plot(lam,I);title('Total Signal');
%FFT
figure(2)
Fs =1/inc;
Fn = Fs/2;
%Total
FTT = fft(I);
FvT = linspace(0, 1, fix(numel(FTT)/2)+1)*Fn;
IvT = 1:numel(FvT);
subplot(3,1,1)
semilogy(FvT, abs(FTT(IvT))*2)
title('Total Signal'); grid; xlim([0 2]);
%First
FT_I1 = fft(I_first);
Fv1 = linspace(0, 1, fix(numel(FT_I1)/2)+1)*Fn;
Iv1 = 1:numel(Fv1);
subplot(3,1,2)
semilogy(Fv1, abs(FT_I1(Iv1))*2)
title('First Signal');grid;xlim([0 2])
%Second
FT_I2 = fft(I_second);
Fv2 = linspace(0, 1, fix(numel(FT_I2)/2)+1)*Fn;
Iv2 = 1:numel(Fv2);
subplot(3,1,3)
semilogy(Fv2, abs(FT_I2(Iv2))*2)
title('Second Signal');grid;xlim([0 2])
% FILTERING
% First Signal
[pk1,loc1]=findpeaks(abs(FT_I1(Iv1)),Fv1);
dn1=0;f1=loc1; %f2=mean(loc1)+dn1;
Wp1 = 1.1*f1/Fn;
Ws1 = 1.01*Wp1;
Rp1 = 1;
Rs1 = 60;
[n11,Wp1] = ellipord(Wp1,Ws1,Rp1,Rs1);
[z,p,k] = ellip(n11,Rp1,Rs1,Wp1);
[sos1,g1] = zp2sos(z,p,k);
figure(3)
freqz(sos1, 2^18, Fs)
set(subplot(2,1,1), 'XLim',Wp1*Fn.*[0.8 1.1])
title('First Signal Magnitude');
set(subplot(2,1,2), 'XLim',Wp1*Fn.*[0.8 1.1])
title('First Signal Phase');
sf1 = smoothdata(filtfilt(sos1, g1, I),'gaussian');
% Second Signal
[pk2,loc2]=findpeaks(abs(FT_I2(Iv2)),Fv2);
dn2=0.0002;f3=mean(loc2)-dn2; f4=mean(loc2)+dn2;
Wp2 = [f3 f4]/Fn;
Ws2 = [0.99 1.01].*Wp2;
Rp2 = 1;
Rs2 = 100;
[n22,Wp2] = ellipord(Wp2,Ws2,Rp2,Rs2);
[z,p,k] = ellip(n22,Rp2,Rs2,Wp2);
[sos2,g2] = zp2sos(z,p,k);
figure(4)
freqz(sos2, 2^18, Fs)
set(subplot(2,1,1), 'XLim',Wp2*Fn.*[0.8 1.1])
title('Second Signal Magnitude');
set(subplot(2,1,2), 'XLim',Wp2*Fn.*[0.8 1.1])
title('Second Signal Phase');
sf2 = smoothdata(filtfilt(sos2, g2, I),'gaussian');
% LENGTH CALCULATION
% First signal length
pp=1;
[p1,l1]=findpeaks(sf1,lam);
FSR1=mean(diff(l1));
FSR11=l1(pp+1)-l1(pp);
length1=0.001*(lam(pp+1)*lam(pp))/(2*n1*FSR11);
%Second signal length
[p2,l2]=findpeaks(sf2,lam);
FSR2=mean(diff(l2));
FSR22=l2(pp+1)-l2(pp);
length2=0.001*(lam(pp+1)*lam(pp))/(2*n2*FSR22);
figure(5)
subplot(4,1,1)
plot(lam,I_first)
title(sprintf(' Original first signal Length = %d',L1/m));grid;
subplot(4,1,2)
plot(lam,sf1)
title(sprintf(' Reconstructed first signal Length = %f',length1));grid;
subplot(4,1,3)
plot(lam,I_second,'r')
title(sprintf(' Original second signal Length = %d',L2/m));grid;
subplot(4,1,4)
plot(lam,sf2,'r')
title(sprintf(' Reconstructed second signal Length = %f',length2));grid;
  1 Comment
Star Strider
Star Strider on 5 Jan 2021
I have lost track of where we are.
  1. I doubt that the passband or stopband ripple magnitudes are significant problems. However also note that the stopband ripple (the first filter is a lowpass filter with a 60 dB stopband rippole ans the second is a bandpass filter with a 100 dB passband ripple) are also the stopband attenuations of the filters. Changing the attenuations would change the signal output with respect to the magnitude of the signal components being eliminated, and the loss of that signal spectral energy would change the filtered output. You would simply need to experiment to determine what works best.
  2. With IIR filters, that is not likely to be a problem, since both of these are significantly shorter than the signal length.
  3. As I mentioned, the frequencies eliminated by the filter will affect the shape of the filtered output. You will need to experiment to get the result you want.
  4. Windowing filters are useful for FIR filters, because the windows correct for a sampled signal having a Fourier transform that does not extend to . I experimented with FIR filters earlier in this discussion, and I encourage you to experiment with them if that is the direction that you want to go. The problem is that they can be very long and therefore ineffecient, so designing them to do what you want and be efficient is an important trade-off. You may have to extend the lengths of your signals to accommodate long filters.
  5. I wish I could. I do not entirely understand what you are doing. The filter design is of course important, and most of my signal processing involves eliminating band-limited noise and unwanted frequencies from physiological signals, rather than reconstructing the signals. Note that eliminating the D-C offset from a signal that does not have a mean of zero (or close to zero) will result in some distortion simply due to the filter algorithms. Normally, that is not a significant problem, however for your signals it appears to be. One option may be to subtract the mean of the signal before you filter it, then add the mean back afterwards. That has worked for me in the past.

Sign in to comment.

Tags

Community Treasure Hunt

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

Start Hunting!