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

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

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

Thank you very much. Do I have to use fftshift even if I use fftfilt function? I am stuck at this point for a couple of days. It would be a great help if you favor me a little bit . I don't know how to use the results of freqz symmterically to each half of the complex ftt result and flip the result.
My pleasure!
Do I have to use fftshift even if I use fftfilt function?
Only if you want to show the result of the filter application. The fftfilt function takes care of everything else.
I don't know how to use the results of freqz symmterically to each half of the complex ftt result and flip the result.’
What I intend by ‘symmetrically’ is to use flipud and then conj on the ‘h’ output of freqz (since the output is a column vector), then apply it to the negative half of the fftshift result. (This is what I meant by this procedure being more trouble than it is worth!) So you would use:
h_lefthalf = flipud(conj(h));
on the left half.
Thank you for your feedback. I think I'm very close but stuck in a loop. I used the fftfilt like the following,
sf=fftfilt(b,I); % b has been calculated using FIR filter.
What is the next step to get the original wavelength domain signal (I_first)?
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=fftfilt(b,I);
Then IFFT to get back I_first?
No. The output of fftfilt (‘sf’ here if I understand your code correctly) is the filtered time-domain signal. The fftfilt function takes care of all those details.
However, If I plot(lam,sf), this graph does not look like I_first? Could you please look at my code one more time and tell me why I'm not getting the same graph?
What I want is:
plot(lam,I_first) and plot(lam,sf) should be same because I used FIR filter only to separate the I_first signal after perfroming FFT. Basically, I'm trying to reconstruct the I_first signal.
It appears that ‘I_First’ has a frequency of about 1/12 cycles/time unit, so if time units are seconds, about 0.08 Hz. So the filter is operating correctly, with respect to that signal, although the filter appears to be adding a constant to it, likely because the signal originally has a mean that is greater than 0. Since lowpass filters (such as this one) are also integrators, and there are other issues with the way it is being filtered, it is integrating the constant offset to produce an increasing amplitude, so all this is entirely consistent.
I base this on a slight variation of the last plot (figure(4) here):
figure(4)
sf=fftfilt(b,I);
plot(lam,sf)
grid
that demonstrates those properties.
So, it appears that your filter is doing what you want, although its design is not optimal with respect to the signal it is filtering.
Thank you so much for your informative explanation. You have been very helpful.
I have added the following lines in the code and got the graphs. Definitely they are not same. You have mentioned the filter design is not optimal. What are the driving parameters that makes the filter design optimal. Why also their amplitude is not same? What are the ways that the reconstructed signal can be very close to the original signal?
figure(4)
sf=fftfilt(b,I); %sf=filter(b,1,I);
subplot(2,1,1)
plot(lam,I_first)
title('Original signal')
subplot(2,1,2)
plot(lam,sf) % wavelength domain signal after filtering
title('Reconstructed signal')
As always, my pleasure!
This is a somewhat difficult signal to filter, however it is possible to do it with an IIR filter and filtfilt, producing a reasonably good result:
Fs = 1/inc;
Fn = Fs/2; % Nyquist Frequency
Wp = 0.1/Fn; % Stopband Frequency (Normalised)
Ws = 1.01.*Wp; % Passband Frequency (Normalised)
Rp = 1; % Passband Ripple
Rs = 20; % Passband Ripple (Attenuation)
[n,Wp] = ellipord(Wp,Ws,Rp,Rs); % Elliptic Order Calculation
[z,p,k] = ellip(n,Rp,Rs,Wp); % Elliptic Filter Design: Zero-Pole-Gain
[sos,g] = zp2sos(z,p,k); % Second-Order Section For Stability
figure(7)
freqz(sos, 2^20, Fs) % Filter Bode Plot
set(subplot(2,1,1), 'XLim',Wp*Fn.*[0.8 1.2]) % Optional
set(subplot(2,1,2), 'XLim',Wp*Fn.*[0.8 1.2]) % Optional
sf2 = filtfilt(sos, g, I); % Filter With IIR Filter
figure(8)
sf=fftfilt(b,I); %sf=filter(b,1,I);
subplot(2,1,1)
plot(lam,I_first)
title('Original signal')
subplot(2,1,2)
plot(lam,sf2) % wavelength domain signal after filtering
title('Reconstructed signal')
I designed a different FIR filter for you to experiment with, and used it with fftfilt. (The only other change in your previous code are the figure numbers.)
This gives a reasonably good result, except for the phase distortion and phase delay:
fcomb = [0.01 0.05 0.14 0.20];
mags = [0 1 0];
dev = [0.5 0.1 0.5];
[n,Wn,beta,ftype] = kaiserord(fcomb,mags,dev,Fs);
hh = fir1(n,Wn,ftype,kaiser(n+1,beta),'noscale');
figure(6)
freqz(hh, 1, 2^20, Fs)
set(subplot(2,1,1), 'XLim',[0 0.5]) % Optional
set(subplot(2,1,2), 'XLim',[0 0.5]) % Optional
figure(7)
sf=fftfilt(hh,I); %sf=filter(b,1,I);
subplot(2,1,1)
plot(lam,I_first)
title('Original signal')
subplot(2,1,2)
plot(lam,sf) % wavelength domain signal after filtering
title('Reconstructed signal')
The phase delay and phase distortion is part of the filtering process, and there is no way to eliminate it using the current functions. (This is the reason filtfilt is preferred, however it is not possible to use it with your signal because of the signal length problem. This filter has 1433 elements. The IIR filter is only 6-order.)
% sf2 = filtfilt(hh, 1, I); % Errors: ‘Data length must be larger than 34371 samples.’
Using filter with this FIR filter produces a reasonably decent result:
sf3 = filter(hh, 1, I);
figure(9)
subplot(2,1,1)
plot(lam,I_first)
title('Original signal')
subplot(2,1,2)
plot(lam,sf3) % wavelength domain signal after filtering
title('Reconstructed signal')
There is some (expected) phase distortion that could likely be prevented by setting the initial conditions of the filter function to something more appropriate than the defaults. Note that I use the FIR filter design procedure that I usually use to design them, letting kaiserord do most of the ‘heavy lifting’. This makes ffor a much easier design process.
.
Awesome. The IIR really works for my signal. Could you please tell me how did you select the following parameters? I can see the resconstructed signal is not much smoooth as like I_first and amplitude is not exactly same.
Wp = 0.1/Fn; % Stopband Frequency (Normalised)
Ws = 1.01.*Wp; % Passband Frequency (Normalised)
Rp = 1; % Passband Ripple
Rs = 20;
freqz(sos, 2^20, Fs)
Most important question:
For example, If I want to reconstruct the I_second signal, what parameters need to be changed here? What would the value of Wp, Ws, Rp and Rs?
What is the problem with the signal length to apply fftfilt? You could vary the lam range as your wish.
Thank you!
First, with respect to the passbands for the filter, I took the Fourier transform of ‘I_first’ and used those results (with a bit of later experimentation) to define the passband. (I originally defined it as a bandpass filter, then thought better of that and designed it as a lowpass filter.)
FT_I1 = fft(I_first);
Fs = 1/inc;
Fn = Fs/2;
Fv = linspace(0, 1, fix(numel(FT_I1)/2)+1)*Fn;
Iv = 1:numel(Fv);
figure(5)
semilogy(Fv, abs(FT_I1(Iv))*2)
grid
xlim([0 1])
To reconstruct ‘I_second’, define the passbands using the same Fourier transform procedure to define its spectral characteristics as I did with ‘I_first’ here. (You will need to change the xlim values appropriately, so comment-out that line until you determine the appropriate frequencies for it and the others.) Then experiment with the appropriate passbands. For it (and the others), a bandpass design will be necessary. Note that this will result in the loss of the D-C (constant) offset, however the shape of the curve itself should be preserved.
If you want to be really adventurous, you can design a stopband-passband FIR filter that will preserve the D-C offset, selectively stop the frequencies you do not want, and selectively pass the frequencies you do want. This is not difficult, however it requires experimentation (therefore time) to get everything correct so that the filter does what you want it to do. You can extend the code I posted to do that.
The problem with the signal length is that I am not certain how you construct the time vector. It would likely not be difficult to extend it, however rather than disrupt the code, I prefer to design filters that work with it as it is. You are certainly welcome to extend it as necessary, then see if the filter I designed will give the correct result with filtfilt and the longer vector.
Thank you very much. I have tried to design the bandpass filter to reconstruct I_second signal. However, the shape of the reconstructed signal is different from the original I_second. I don't know why the FFT of I_second gives two peaks although it is a single cosine signal. The FFT of I_third is similar to I_second but slightly shifted. Could you please tell me how to select RP and Rs? I can see Rs has a huge impact on the reconstructed signal. Please tell me where I made a mistake in the follwoing code that is not providing me the reconstructed signal similar to original signal.
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
figure(2)
Fs = 1/inc;
Fn = Fs/2; % Nyquist Frequency
FT_I1 = fft(I_second);
Fs =1/inc;
Fn = Fs/2;
Fv = linspace(0, 1, fix(numel(FT_I1)/2)+1)*Fn;
Iv = 1:numel(Fv);
semilogy(Fv, abs(FT_I1(Iv))*2)
grid
xlim([0 2])
% Filtering
Wp = [1.17 1.22]/Fn; % Stopband Frequency (Normalised)
Ws = [0.98 1.02].*Wp; % Passband Frequency (Normalised)
Rp = 1; % Passband Ripple
Rs = 20; % Passband Ripple (Attenuation)
[n,Wp] = ellipord(Wp,Ws,Rp,Rs); % Elliptic Order Calculation
[z,p,k] = ellip(n,Rp,Rs,Wp); % Elliptic Filter Design: Zero-Pole-Gain
[sos,g] = zp2sos(z,p,k); % Second-Order Section For Stability
figure(3)
freqz(sos, 2^20, Fs) % Filter Bode Plot
set(subplot(2,1,1), 'XLim',Wp*Fn.*[0.8 1.1]) % Optional
set(subplot(2,1,2), 'XLim',Wp*Fn.*[0.8 1.1]) % Optional
sf2 = filtfilt(sos, g, I); % Filter With IIR Filter
figure(4)
subplot(2,1,1)
plot(lam,I_second)
title('Original signal')
subplot(2,1,2)
plot(lam,sf2) % wavelength domain signal after filtering
title('Reconstructed signal')
When you combined (added) all the ‘I’ signals together, you ended up with a modulated carrier, specifically a double-siideband suppressed-carrier signal. (The Wikipedia article discussed amplitude modulation in much more detail.) It is a bit complicated to go into here, however that explains the two peaks and the reconstructed signal. You get much closer to the original signal (although with a significantly reduced amplitude) with:
Wp = [1.216 1.217]/Fn;
Amplitude modulation also introduces a small amount of frequency (or phase) modulation, the reason the phase is not exactly the same as that of the original signal.
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)
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)

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)
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

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!