Why isn't low pass filter centered around zero?
18 views (last 30 days)
Show older comments
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
0 Comments
Accepted Answer
Paul
on 21 Dec 2022
Edited: Paul
on 14 Jun 2023
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 Comments
Paul
on 26 Dec 2022
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.
Paul
on 27 Dec 2022
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.
More Answers (0)
See Also
Categories
Find more on Fourier Analysis and Filtering in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!