Why isn't low pass filter centered around zero?

12 views (last 30 days)
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.
x1 = -8.4645e+03
x2 = 8.2661e+03
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.
x1 = -1.0365e+06
x2 = 6.2328e+05
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

Accepted Answer

Paul
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);
Nfft = 121
myplot(xFreq, NUMfft); title('Fs = 48 KHz') % slightly off-centered; x1 should be -x2.
x1 = -8.3306e+03
x2 = 8.3306e+03
Fs = 100e6; % my number
Fp = 500e3; % my number
[NUMfft, xFreq] = getFilter(N, Fp, Fs, Rp, Rst);
Nfft = 121
myplot(xFreq, NUMfft); title('Fs = 100 MHz') % very off-centered; x1 should be -x2.
x1 = -8.2645e+05
x2 = 8.2645e+05
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
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))
xval = 1×7
1 2 3 4 5 6 7
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))
X_c(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.

Sign in to comment.

More Answers (0)

Products


Release

R2020a

Community Treasure Hunt

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

Start Hunting!