FFT-filter-IFFT problem
Nicolas
I'm studying the effects of a bandpass filter on a single rectangular pulse. I've measured how the transference changes with the frequency (function 'Transfer' in the code) and also the phase difference between the input and output signals with frequency (function 'Phase' in the code).
The code is as follows:
clear all
% Characterization of the input signal
fs = 1e12;
ts = 1/fs;
T=1e-8;
N=2*T*fs;
t = -T:ts:T;
tau=T/10;
Pulse = rectpuls(t,tau);
% FFT of the input
M=2^nextpow2(N);
f=linspace(-fs/2,fs/2,M);
x=fft(Pulse,M);
sizef=size(f,2);
absx=abs(x(1:sizef));
anglex=angle(x(1:sizef));
% Effects of the bandpass circuit
Transfer=zeros(1,sizef);
Phase=zeros(1,sizef);
for i=1:sizef
Transfer(i)=(sqrt(3.24*10^(-7)*(f(i))^2+(1.84*10^4)*(f(i))^(-2)+3.4))^(-1);
Phase(i)=atan(f(i)*2.93*10^(-4)-56.16/f(i));
end
absy=absx.*Transfer;
angley=anglex+Phase;
% IFFT to get back to time domain
y=absy.*cos(angley)+1j*absy.*sin(angley);
iffty=ifft(y,M);
subplot(2,1,1)
plot(t,Pulse,'-b');
axis([-T T -1 2]);
title('Input pulse');
xlabel('Time [seg]');
ylabel('Voltage [V]');
subplot(2,1,2)
plot(t,real(iffty(1:size(t,2))),'-r')
title('Output pulse')
xlabel('Time [seg]');
ylabel('Voltage [V]');
*******************************************
As you can see, there is a strange effect in the central zone of the plot of the output pulse (the region of interest). It's like the code is plotting twice and then just join both group of points. I'd like any idea to solve this issue, because the plot can't be like that.
If anyone has any suggestion about the code, i'll be appreciated to receive it.
dbd
> Hello.
> ...
> If anyone has any suggestion about the code, i'll be appreciated to receive it.
First, you need to look at a basic description of the discrete digital
process you are trying to accomplish. Start with something like:
http://www.dspguide.com/ch18/2.htm
(Try the rest of the book as well. You can get it on the web. Free.)
Second, the code looks like an attempt to obfuscate by writing bad
Fortran in Matlab. It probably isn't important to use a power of 2
transform size (but it is necessary to use a transform size larger
than the impulse length of the signal plus the impulse length of the
filter function minus 1; if this doesn't make sense, read more of the
referenced book or some other resources). It is potentially confusing
in Matlab to redefine "i" from the default value as the square root of
-1; particularly in the context of phase calculations.
Third, The fft function calculates the discrete Fourier transform
(DFT). It has some similar, but not identical properties to the
(continuous/infinite) Fourier transform. Blindly mapping the
properties of the later to the former is commonly attempted by Matlab
beginners with rare success. Outside of symbolic manipulation, you are
limited to finite discrete processes and their characteristics and
limitations.
Dale B. Dalrymple
Greg Heath
> Hello.
>
> I'm studying the effects of a bandpass filter on a single rectangular pulse. I've measured how the transference changes with the frequency (function 'Transfer' in the code) and also the phase difference between the input and output signals with frequency (function 'Phase' in the code).
>
> The code is as follows:
> As you can see, there is a strange effect in the central zone of the plot of the output pulse (the region of interest). It's like the code is plotting twice and then just join both group of points. I'd like any idea to solve this issue, because the plot can't be like that.
> If anyone has any suggestion about the code, i'll be appreciated to receive it.
See my recent post
Filtering using MATLAB's ifftshift
http://groups.google.com/group/comp.soft-sys.matlab/msg/ac495478ad970218?hl=en
Hope this helps.
Greg
Nicolas
Thanks again.
dbd
> On May 14, 8:36 pm, "Nicolas " <nico_spect...@hotmail.com> wrote:> Hello.
> > If anyone has any suggestion about the code, i'll be appreciated to receive it.
>
> See my recent post
>
> Filtering using MATLAB's ifftshift
>
>
> Hope this helps.
>
> Greg
Unfortunately, this interesting post on Matlab syntax fails to discuss
any of the requirements for correctly performing filtering via fft and
ifft. Correct calculation does not require "shift" syntax. (When you
understand fft-filter-ifft, go to whatever syntax you prefer.)
If you are interested in the filtering process and what it does
require, regardless of the syntax of the implementation, here is a
simplified example of such a process in Matlab presented in a format
consistent with the book I referenced:
http://www.dspguide.com/ch18/2.htm
% Section 1
% An fft-filter-ifft example, pulse and lowpass filter
clear all
% The filter: lowpass, a truncated sinc
M = 14; % an even number for convenience next below
filt = sinc(0.2*(-M/2:-1+M/2));
% The signal: a single pulse
N = 20; % an even number for convenience next below
signal = 1.5*[ones(1,N/2) zeros(1,N/2)];
% DFT the filter kernel and signal
% The transform size must be at least M + N -1
fftfilt = fft(filt,M+N);
fftsignal = fft(signal,M+N);
% The Inverse DFT of the frequency domain filter process
filtsignal = real(ifft(fftsignal .* fftfilt));
figure
plot(1:M+N,filtsignal,'r*', ...
1:M,filt,'ko',1:N,signal,'bx');grid on;
axis([0 M+N+1 -1.15 9.5])
title('Filter a FINITE discrete sequence with a discrete filter')
xlabel('Finite Data sequence: x Filter sequence: o Output: *')
% Section 2
% For the case of an infinite periodic data sequence
% The following code works for N >= M
fftfilt = fft(filt,N);
fftsignal = fft(signal,N);
% The Inverse DFT of the frequency domain filter process
filtsignal = real(ifft(fftsignal .* fftfilt));
figure
plot(1:N,filtsignal,'r*', ...
1:M,filt(1:M),'ko',1:N,signal,'bx');grid on;
axis([0 N+1 -1.15 9.5])
title('Filter a PERIODIC discrete sequence with a discrete filter')
xlabel('One Period Data sequence: x Filter sequence: o One Period
of Output: *')
The first section is the basic convolution operation via fft-ifft. The
FFT size must be large enough to span the sizes of the filter and data
sequences minus 1. This is true of time domain convolution as well.
The book chapter cited gives the details of applying this process to
long data sequences. The additional processing for long data sequences
takes forms such as "overlap-add".
One exception to the FFT size requirement is in the case of a known
periodic data sequence. This is illustrated in the second section of
the code. In the periodic sequence case, the aliasing due to circular
convolution from using too small a transform for linear convolution
actually performs the "overlap-add" operations. Note that the results
of the second section are not a subset of the results of the first
section.
Good luck
Dale B. Dalrymple