# Why isn't low pass filter centered around zero?

11 ビュー (過去 30 日間)
Paul Hoffrichter 2022 年 12 月 20 日

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

サインインしてコメントする。

### 採用された回答

Paul 2022 年 12 月 21 日

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 件のコメント7 件の古いコメントを表示7 件の古いコメントを非表示
Paul 2022 年 12 月 26 日
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))')
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))')
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))')
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 2022 年 12 月 27 日
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')
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.

サインインしてコメントする。

### カテゴリ

Help Center および File ExchangeCalculus についてさらに検索

R2020a

### Community Treasure Hunt

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

Start Hunting!

Translated by