Why isn't low pass filter centered around zero?

14 vues (au cours des 30 derniers jours)
Paul Hoffrichter
Paul Hoffrichter le 20 Déc 2022
Modifié(e) : Paul le 26 Juin 2023
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

Réponse acceptée

Paul
Paul le 21 Déc 2022
Modifié(e) : Paul le 14 Juin 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 commentaires
Paul
Paul le 26 Déc 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).
Now, let's compute the 7-point DFT of x[n] via fft and add it to the plot.
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
Paul le 27 Déc 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.

Connectez-vous pour commenter.

Plus de réponses (0)

Produits


Version

R2020a

Community Treasure Hunt

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

Start Hunting!

Translated by