Why isn't low pass filter centered around zero?
Afficher commentaires plus anciens
I tried several MATLAB filter examples (in R2020a), but all of them appear to be off-centered. Below is example taken from https://www.mathworks.com/help/dsp/ug/lowpass-filter-design.html
I used this example to try to create a low pass filter for my application. But the 2-sided fft plot of the real filter numerator appears to be off-centered especially when I replace the example numbers with my own numbers. Why isn't the fft magnitude symmetric around x=0?
N=120;
Ap = 0.01;
Ast = 80;
Rp = (10^(Ap/20) - 1)/(10^(Ap/20) + 1);
Rst = 10^(-Ast/20);
Fs = 48e3; % from example
Fp = 8e3; % from example
[NUMfft, xFreq] = getFilter(N, Fp, Fs, Rp, Rst);
myplot(xFreq, NUMfft); title('Fs = 48 KHz') % slightly off-centered; x1 should be -x2.
Fs = 100e6; % my number
Fp = 500e3; % my number
[NUMfft, xFreq] = getFilter(N, Fp, Fs, Rp, Rst);
myplot(xFreq, NUMfft); title('Fs = 100 MHz') % very off-centered; x1 should be -x2.
function [NUMfft, xFreq] = getFilter(N, Fp, Fs, Rp, Rst)
NUM = firceqrip(N,Fp/(Fs/2),[Rp Rst],'passedge');
NUMfft = abs(fftshift(fft(NUM)));
xFreq = linspace( -1, 1-1/length(NUMfft), length(NUMfft) ) * Fs/2;
end
function myplot(x, y)
thresh = 0.8;
plot(x,y, '-b.')
xline(0,'r')
v1 = find(y>thresh, 1, 'first');
v2 = find(y>thresh, 1, 'last');
x1 = x(v1)
x2 = x(v2)
xline(x1,'k'), xline(x2,'k')
ylim([0.8 1])
end
Réponse acceptée
Hi Paul,
The fft length is odd, so the computation of xFreq needs to be as shown below.
N=120;
Ap = 0.01;
Ast = 80;
Rp = (10^(Ap/20) - 1)/(10^(Ap/20) + 1);
Rst = 10^(-Ast/20);
Fs = 48e3; % from example
Fp = 8e3; % from example
[NUMfft, xFreq] = getFilter(N, Fp, Fs, Rp, Rst);
myplot(xFreq, NUMfft); title('Fs = 48 KHz') % slightly off-centered; x1 should be -x2.
Fs = 100e6; % my number
Fp = 500e3; % my number
[NUMfft, xFreq] = getFilter(N, Fp, Fs, Rp, Rst);
myplot(xFreq, NUMfft); title('Fs = 100 MHz') % very off-centered; x1 should be -x2.
function [NUMfft, xFreq] = getFilter(N, Fp, Fs, Rp, Rst)
NUM = firceqrip(N,Fp/(Fs/2),[Rp Rst],'passedge');
NUMfft = abs(fftshift(fft(NUM)));
% xFreq = linspace( -1, 1-1/length(NUMfft), length(NUMfft) ) * Fs/2;
Nfft = numel(NUMfft)
xFreq = ((-(Nfft-1)/2) : (Nfft-1)/2)/Nfft * Fs;
end
function myplot(x, y)
thresh = 0.8;
plot(x,y, '-b.')
xline(0,'r')
v1 = find(y>thresh, 1, 'first');
v2 = find(y>thresh, 1, 'last');
x1 = x(v1)
x2 = x(v2)
xline(x1,'k'), xline(x2,'k')
ylim([0.8 1.1])
end
9 commentaires
Paul Hoffrichter
le 22 Déc 2022
Paul
le 22 Déc 2022
HTH. Also, if the fft length is even, then we'd have
xFreq = (-Nfft/2 : (Nfft/2-1))/Nfft * Fs;
Paul Hoffrichter
le 22 Déc 2022
Here's an explanation that justifies the expressions for xFreq. It's not the only way to explain it, and other ways might be better in relation to the continuous time Fourier transform, but hopefully this approach gives an intuitive explanation based on what the output of fft is and how it relates to the ifft.
Let's define a finite duration discrete-time signal of length N = 7
N = 7;
x = @(n) (n+1).*(n >= 0 & n <= N-1);
Let's plot the signal for a few values of n. By definition, the sequence is is zero for n < 0 and n >= N.
n = -3:10;
figure
stem(n,x(n))
xlabel('n');
ylabel('x[n]')
axis('padded')
Let's compute the Discrete Fourier Tranform (DFT) of the finite duration portion of the signal that we care about, i.e., for 0 <= n <= N-1. The function mydft computes the DFT based on its definiton and is instructive, though probably not the most efficient nor probably the most accurate.
n = 0:(N-1);
Xdft = mydft(x(n));
Let's compare Xdft to the output of fft, which also computes the DFT
norm(Xdft - fft(x(n)))
So it looks like mydft correctly computes the DFT.
Now, let's look a bit more carefully at how the DFT is computed. The function mydft shows that each element of Xdft is computed from a corresponding angle in the vector theta. Each element of theta is used in a complex exponential and therefore represents an angle around the unit circle.
Rewrite theta here using sym so that we can better see each value.
theta = sym(0:(N-1))/sym(N)*sym(2)*sym(pi)
The angles in theta that correspond to the elements of Xdft are equally spaced around the unit circle with dtheta = 2*pi/N. The last point is NOT 2*pi, it is one dtheta short of 2*pi.
Let's make a table that shows the corespondence between theta and Xdft
table(theta.' , Xdft.','VariableNames',{'theta' , 'Xdft'})
Now, let's consider going in the other direction, i.e., recovering the original sequence from its DFT. The function myidft computes the Inverse DFT (IDFT) based on its definition. It recovers the original N-point sequence.
xval = myidft(Xdft);
xval.'
Inspection of myidft shows that each element of xval is computed based on a corresponding angle in the vector theta, where the angles in theta in the IDFT are the exact same angles used to compute the DFT.
Furthermore, in the reconstruction of each term in the IDFT, the angle in theta is multiplied by integers in the complex exponential. Therefore, we can take any angle in theta and add or subtract an integer multiple of 2*pi and have no effect on the IDFT computation. Let's subtract 2*pi from each element in theta that is greater than or equal to pi
theta(theta >= sym(pi)) = theta(theta >= sym(pi)) - 2*sym(pi)
This new theta vector still corresponds to the same elemnts in Xdft. Let's update our table
table(theta.', Xdft.','VariableNames',{'theta' , 'Xdft'})
The table looks odd because theta is no longer monotonic. Let's fix that by sorting theta and then rearranging Xdft so that each element of Xdft maintains the correct correspondence to its element of theta after sorting
[~,index] = sort(double(theta));
theta = theta(index)
Xdft = Xdft(index);
theta is now symmetric around zero.
Finally, let's apply fftshift to the output of fft and create our table
Xfftshift = fftshift(fft(x(n)));
table(theta.', Xdft.', Xfftshift.', 'VariableNames',{'theta' , 'Xdft', 'Xfftshift'})
It should be clear that our current theta vector corresponds to the elements of the output of fftshift, because they are the same as as the elements of our (rearranged) Xdft
This vector theta can also be computed as
sym((-(N-1)/2) : ((N-1)/2))/sym(N)*sym(2)*sym(pi)
Recall that we started with N = 7, which is an odd number. How would things change if N is even? Let's try N = 8, and compute the theta vector corresponding to the DFT
N = 8;
theta = sym(0:(N-1))/sym(N)*sym(2)*sym(pi)
Note that pi shows up exactly in the middle of the array. If we follow the development above by subtracting 2*pi from each element greater than or equal to pi and sorting we get
theta(theta >= sym(pi)) = theta(theta >= sym(pi)) - 2*sym(pi);
[~,index] = sort(double(theta));
theta = theta(index)
We see that -pi shows up as the first element on the left. This result is consistent with the convention used by fftshift when N is even. This vector theta can be computed directly as
sym((-N/2) : (N/2 - 1))/sym(N)*sym(2)*sym(pi)
In summary we have the angles that correpond to the output of fftshift are given by
For N odd:
theta = ((-(N-1)/2) : ((N-1)/2))/N*2*pi (rad, or rad/sample)
For N even:
theta = ((-N/2) : (N/2 - 1))/N*2*pi (rad, or rad/sample)
If we don't include the 2*pi factor in theta, the units effectively change to cycles, or cycles/sample, keeping in mind that the 2*pi factor would still need to be included in the complex exponential in the DFT and IDFT computations (because the angle input to exp always has to be in radians).
The last aspect of the problem is to consider the effect of the sampling frequency if the elements of x represent equally spaced samples of a continuous-time signal. Let Ts be the sampling period and Fs = 1/Ts be the sampling frequency.
In short, we can define a vector of frequencies for the DFT and IDFT from theta, without the 2*pi subtraction, as
omega = theta/Ts (rad/sec)
With this definiton, the DFT would be computed with the argument in the complex exponential as omega*Ts. Of course, the end result for Xdft would be the same, but we can view each element of Xdft as corresponding to an element of omega. However, we have to keep in mind that each element of omega still corresponds to an angle via theta = omega*Ts. So we can still subtract the 2*pi/Ts from selected elements of omega. In short, the end result would be:
For N odd:
omega = ((-(N-1)/2) : ((N-1)/2))/N*2*pi/Ts (rad/sec)
For N even:
omega = ((-N/2) : (N/2 - 1))/N*2*pi/Ts (rad/sec)
Note that with N even, the left most element in omega will be -pi/Ts, which is the negative of the Nyquist frequency.
By the same argument we can replace 2*pi/Ts with Fs in the definition of omega, in which case omega would have units of Hz, and then the DFT and IDFT would include a factor of 2*pi*Ts in the complex exponential.
function Z = mydft(z)
N = numel(z);
Z = zeros(1,N);
theta = (0:(N-1))/N*2*pi;
n = (0:(N-1));
for kk = 0:(N-1)
Z(kk+1) = sum(z.*exp(-1j*theta(kk+1).*n));
end
end
function z = myidft(Z)
K = numel(Z);
z = zeros(1,K);
theta = (0:(K-1))/K*2*pi;
k = (0:(K-1));
for nn = 0:(K-1)
z(nn+1) = sum(Z.*exp(1j*theta(nn+1).*k))/K;
end
end
Nice well written tutorial. I stepped through your code in a .mlx file, and I figured out the equivalent linspace versions. Thanks again.
N = 7;
sym( -(N-1)/2 : (N-1)/2 )/N * 2*pi % N is odd
sym( linspace( -1+1/N, 1-1/N, N)*pi ) % N is odd
N = 8;
sym( ( -N/2 : (N/2-1) )/N * 2*pi ) % N is even
sym( linspace( -1, 1-2/N, N)*pi ) % N is even
A general expression that works for N even or odd is
ceil(linspace(-N/2,N/2-1,N))/N*2*pi
N = 7;
sym(ceil(linspace(-N/2,N/2-1,N))) /N*2*pi
N = 8;
sym(ceil(linspace(-N/2,N/2-1,N))) /N*2*pi
Excellent! Turns the 5-line if-else-end into a single line when need to handle both even and odd! Thanks again! To summarize for handy future reference (especially when code-reviewing):
N = 7; % N odd
th1 = sym( ceil( linspace(-N/2, N/2-1, N) ) )/N * 2*pi % odd or even
th2 = sym ( ceil( -N/2 : N/2-1 ) )/N * 2*pi % odd or even
th3 = sym( -(N-1)/2 : (N-1)/2 )/N * 2*pi % odd only
th4 = sym( linspace( -1+1/N, 1-1/N, N) ) * pi % odd only
N = 8 % N even
phi1 = sym( ceil( linspace(-N/2, N/2-1, N) ) )/N * 2*pi % odd or even
phi2 = sym( ceil( -N/2 : N/2-1 ) )/N * 2*pi % odd or even
phi3 = sym( ( -N/2 : (N/2-1) ) )/N * 2*pi % N is even
phi4 = sym( linspace( -1, 1-2/N, N) ) * pi % N is even
% Check - Odd PASSED:
format long
N = 121;
th1 = ( ceil( linspace(-N/2, N/2-1, N) ) )/N * 2*pi; % odd or even
th2 = ( ceil( -N/2 : N/2-1 ) )/N * 2*pi; % odd or even
th3 = ( -(N-1)/2 : (N-1)/2 )/N * 2*pi; % odd only
th4 = ( linspace( -1+1/N, 1-1/N, N) ) * pi; % odd only
[ max(th1-th2) max(th1-th3) max(th1-th4) ]
plot(th1-th4)
% Check - Even PASSED:
N=120
phi1 = ( ceil( linspace(-N/2, N/2-1, N) ) )/N * 2*pi; % odd or even
phi2 = ( ceil( -N/2 : N/2-1 ) )/N * 2*pi; % odd or even
phi3 = ( ( -N/2 : (N/2-1) ) )/N * 2*pi; % N is even
phi4 = ( linspace( -1, 1-2/N, N) ) * pi; % N is even
[ max(phi1-phi2) max(phi1-phi3) max(phi1-phi4) ]
plot(phi1-phi4)
Happy New Year!
Happy New Year to you!
Here's a bit more information that might be of interest that may give a little more insight into what fftshift is doing and why our expession for the frequency vector corresponding to the the post-fftshifted result is what we want.
Again, define our length-7 signal x[n]:
N = 7;
x = @(n) (n+1).*(n >= 0 & n <= N-1);
Now we compute the Discrete Time Fourier Transform (DTFT) of x[n]. The DTFT is a function of the discrete time angular frequency w (rad, or rad/sample if you prefer), which itself covers the entire real line. The DTFT is periodic, and its period is 2pi. Because x[n] = 0 for n < 0, and it is finite duration (x[n] = 0 for n >= N), we can use freqz to compute its DTFT. Here, we compute the DTFT (X(w)) for a couple of periods, from -3pi to +3pi
wdtft = linspace(-3*pi,3*pi,10000);
xval = x(0:(N-1));
Xdtft = freqz(xval,1,wdtft);
figure
haxmag = subplot(211);hold on;plot(haxmag,wdtft,abs(Xdtft));ylabel('abs(X(\omega))')
haxphs = subplot(212);hold on;plot(haxphs,wdtft,angle(Xdtft));ylabel('angle(X(\omega))')
xlabel('\omega rad/sample')
All of the information about x[n] is encoded in a single period of X(w), and x[n] can be recovered from any interval of X(w) that covers 2pi (discussed in the DTFT link above).
Xdft = fft(xval);
wdft = (0:(N-1))/N*2*pi;
stem(haxmag,wdft,abs(Xdft));
stem(haxphs,wdft,angle(Xdft));
The elements of the DFT are frequency domain samples of the DTFT. This result is not a coincidence. Furthermore, the DFT elements using the convention of fft live in the one-period interval of the DTFT for 0 <= w <= 2pi. That's fine, because all of the information about x[n] is encoded in any 2pi interval of the DTFT.
However, it may be more appealing to consider the interval -pi <= w <= pi so that we get a symmetric look to the plot. In this case, we use fftshift to modify the output of fft so that the DFT elements corresponding to w >= pi are left-right flipped and shifted to the left side of the output, and we redefine our frequency vector as discussed above.
figure
haxmag = subplot(211);hold on;plot(haxmag,wdtft,abs(Xdtft));ylabel('abs(X(\omega))')
haxphs = subplot(212);hold on;plot(haxphs,wdtft,angle(Xdtft));ylabel('angle(X(\omega))')
xlabel('\omega rad/sample')
Xdft = fftshift(fft(xval));
wdft = ceil(-(N/2):(N/2-1))/N*2*pi;
stem(haxmag,wdft,abs(Xdft));
stem(haxphs,wdft,angle(Xdft));
Again, the output from fftshift are samples of the DTFT, as must be the case. Of course, we haven't changed anything on the plot for w >=0, and the frequency spacing between the DFT samples hasn't changed either.
Now, let's do the same thing with a 16-point DFT of our 7-point sequence. We haven't changed the underlying signal x[n], so there's on need to recompute its DTFT.
figure
haxmag = subplot(211);hold on;plot(haxmag,wdtft,abs(Xdtft));ylabel('abs(X(\omega))')
haxphs = subplot(212);hold on;plot(haxphs,wdtft,angle(Xdtft));ylabel('angle(X(\omega))')
xlabel('\omega rad/sample')
Nfft = 16;
Xdft = fftshift(fft(xval,Nfft));
wdft = ceil(-(Nfft/2):(Nfft/2-1))/Nfft*2*pi;
stem(haxmag,wdft,abs(Xdft));
stem(haxphs,wdft,angle(Xdft));
These plots show that zero-padding the FFT yields finer spacing of the frequency domain sampling of the DTFT. The plot also shows that our calculation of wdft for an even-length FFT after fftshift is still correct because the fftshfited DFT elements are samples of the DTFT at the frequencies in wdft.
Let's complete the story by adding the continuous-time Fourier transform (CTFT) and sampling into the analysis.
Define a continuous-time signal
syms t w real
x_c(t) = (1+t*4)*rectangularPulse(-0.125,1.625,t);
Plot the signal
hx = gca;hold on
fplot(hx,x_c,[-1 10])
Assuming a sampling period of Ts = 0.25
Ts = 0.25;
Take 7 samples starting from t = 0.
N = 7;
n = (0:(N-1));
xval = double(x_c(n*Ts))
These are the same samples of x[n] that we've been working with. Also, x_c(nT) is zero for n < 0 and n > 6. Add the samples to the time domain plot
stem(hx,n*Ts,xval)
Compute the CTFT of x_c(t)
X_c(w) = simplify(fourier(x_c(t),t,w))
Plot the CTFT. It looks a lot like the interval of the DTFT of x[n] between -pi and pi
figure
haxmag = subplot(211);hold on;fplot(abs(X_c(w)),[-20 20]),grid,ylabel('abs')
haxphs = subplot(212);hold on;fplot(angle(X_c(w)),[-20 20]),grid,ylabel('angle')
xlabel('\omega (rad/sec)')
Compute the interval of the DTFT of x[n] from -pi to pi
wdtft = linspace(-pi,pi,1000);
Xdtft = freqz(xval,1,wdtft);
Add the DTFT to the plot. Scale the frequency vector by 1/Ts to convert to rad/sec and scale the DTFT by Ts to match the CTFT.
plot(haxmag,wdtft/Ts,abs(Xdtft)*Ts);
plot(haxphs,wdtft/Ts,angle(Xdtft));
Now compute the DFT, fftshift, and add to the plot. Use the formula to compute the fftshifted DFT frequency samples in rad and scale by 1/Ts to convert to rad/sec. Scale the DFT by Ts to match the scaled DTFT.
Xdft = fftshift(fft(xval));
wdft = ceil(-(N/2):(N/2-1))/N*2*pi/Ts;
stem(haxmag,wdft,abs(Xdft)*Ts,'k');
stem(haxphs,wdft,angle(Xdft),'k');
DTFT*Ts approximates the CTFT for -wn < w < wn, where wn is the Nyquist frequency, and Ts*DFT are frequency domain samples of the sclaed DTFT over that same interval, so the fftshifted DFT*Ts are approximate frequency domain samples of the CTFT.
Plus de réponses (0)
Catégories
En savoir plus sur Parametric Spectral Estimation dans Centre d'aide et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
