This example showcases zoom FFT, which is a signal processing technique used to analyze a portion of a spectrum at high resolution. DSP System Toolbox™ offers this functionality in MATLAB™ through the
dsp.ZoomFFT System object, and in Simulink through the zoom FFT library block.
A signal's resolution is bounded by its length. We will illustrate this fact by a simple example. Consider a signal formed by the sum of two sine waves:
L = 32; % Frame size Fs = 128; % Sample rate res = Fs/L; % Frequency resolution f1 = 40; % First sine wave frequency f2 = f1 + res; % Second sine wave frequency sn1 = dsp.SineWave('Frequency',f1,'SampleRate',Fs,'SamplesPerFrame',L); sn2 = dsp.SineWave('Frequency',f2,'SampleRate',Fs,'SamplesPerFrame',L,... 'Amplitude',2); x = sn1() + sn2();
Let's compute the FFT of x and plot the magnitude of the FFT. Note that the two sine waves are in adjacent samples of the FFT. This means that you cannot discriminate between frequencies closer than Fs/L.
X = fft(x); stem(Fs/L*(0:length(X)-1)-Fs/2,abs(fftshift(X))/L) grid on; xlabel('Frequency (Hz)') ylabel('Magnitude') title('Two-sided spectrum')
Suppose you have an application for which you are only interested in a sub-band of the Nyquist interval. The idea behind zoom FFT is to retain the same resolution you would achieve with a full-size FFT on your original signal by computing a small FFT on a shorter signal. The shorter signal comes from decimating the original signal. The savings come from being able to compute a much shorter FFT while achieving the same resolution. This is intuitive: for a decimation factor of D, the new sampling rate is Fsd = Fs/D, and the new frame size (and FFT length) is Ld = L/D, so the resolution of the decimated signal is Fsd/Ld = Fs/L.
DSP System Toolbox offers zoom FFT functionality for MATLAB and Simulink, through the dsp.ZoomFFT System object and the zoom FFT library block, respectively. The next sections will discuss the zoom FFT algorithm in more detail.
Before discussing the algorithm used in dsp.FFT, we present the mixer approach, which is a popular zoom FFT method.
For the example here, assume you are only interested in the interval [16 Hz, 48 Hz].
BWOfInterest = 48 - 16; Fc = (16 + 48)/2; % Center frequency of bandwidth of interest
The achievable decimation factor is:
BWFactor = floor(Fs/BWOfInterest);
The mixer approach consists of first shifting the band of interest down to DC using a mixer, and then performing lowpass filtering and decimation by a factor of BWFactor (using an efficient polyphase FIR decimation structure).
Let's design the decimation filter's coefficients using the function
B = designMultirateFIR(1,BWFactor); D = dsp.FIRDecimator('DecimationFactor',BWFactor,'Numerator',B);
Now, let's mix the signal down to DC, and filter it through the FIR decimator:
for k = 1:10 % Run a few times to eliminate transient in filter x = sn1()+sn2(); % Mix down to DC indVect = (0:numel(x)-1).' + (k-1) * size(x,1); y = x .* exp(-2*pi*indVect*Fc*1j/Fs); % Filter through FIR decimator xd = D(y); end
Let's now take the FFT of the filtered signal (note that the FFT length is reduced by BWFactor, or the decimation length, compared to regular FFT, while maintaining the same resolution):
fftlen = numel(xd); Xd = fft(xd); figure Ld = L/BWFactor; Fsd = Fs/BWFactor; F = Fc + Fsd/fftlen*(0:fftlen-1)-Fsd/2; stem(F,abs(fftshift(Xd))/Ld) grid on xlabel('Frequency (Hz)') ylabel('Magnitude') title('Zoom FFT Spectrum. Mixer Approach.')
The complex-valued mixer adds an extra multiplication for each high-rate sample, which is not efficient. the next section presents an alternative, more efficient, zoom FFT approach.
An alternative zoom FFT method takes advantage of a known result from bandpass filtering (also sometimes called under-sampling): Assume we are interested in the band [F1,F2] of a signal with sampling rate Fs Hz. If we pass the signal through a complex (one-sided) bandpass filter centered at Fc = (F1+F2)/2 and with bandwidth BW = F2 - F1, and then downsample it by a factor of D = floor(Fs/BW), we will bring down the desired band to baseband.
In general, if Fc cannot be expressed in the form k*Fs/D (where K is an integer), then the shifted, decimated spectrum will not be centered at DC. In fact, the center frequency Fc will be translated to :
Fd = Fc - (Fs/D)*floor((D*Fc + Fs/2)/Fs)
In this case, we can use a mixer (running at the low sample rate of the decimated signal) to center the desired band to zero Hertz.
Using the example from the previous section, we obtain the coefficients of the complex bandpass filter by modulating the coefficients of the designed lowpass filter:
Bbp = B .*exp(1j*(Fc / Fs)*2*pi*(0:numel(B)-1)); D = dsp.FIRDecimator('DecimationFactor',BWFactor,'Numerator',Bbp);
Let's now perform the filtering and the FFT:
for k = 1:10 % Run a few times to eliminate transient in filter x = sn1()+sn2(); xd = D(x); end Xd = fft(xd); figure stem(F,abs(fftshift(Xd))/Ld) grid on xlabel('Frequency (Hz)') ylabel('Magnitude') title('Zoom FFT Spectrum. Bandpass Sampling Approach.')
The filter from the previous section is a single-stage filter. We can achieve a less computationally complex filter by using a multistage design instead. This is the approach followed in dsp.ZoomFFT.
Let's consider the following example, where the input sample rate is 48 kHz, and the bandwidth of interest is the interval [1500,2500] Hz. The achievable decimation factor is then floor(48000/1000) = 48.
Let's first design a single-stage decimator:
Fs = 48e3; Fc = 2000; % Bandpass filter center frequency BW = 1e3; % Bandwidth of interest Ast = 80; % Stopband attenuation BWFactor = floor(Fs/BW); B = designMultirateFIR(1,BWFactor,12,80); Bbp = B .*exp(1j*(Fc / Fs)*2*pi*(0:numel(B)-1)); D_single_stage = dsp.FIRDecimator('DecimationFactor',BWFactor,'Numerator',Bbp);
Now, let's design the filter using a multistage design, while maintaining the same stopband attenuation and transition bandwidth as the single-stage case (see
kaiserord for details on the transition width computation):
tw = (Ast - 7.95) / ( numel(B) * 2.285); fmd = fdesign.decimator(BWFactor,'Nyquist',BWFactor,... 'TW,Ast', tw * Fs / (2*pi),Ast,Fs); D_multi_stage = multistage(fmd,'HalfbandDesignMethod','equiripple','SystemObject',true); fprintf('Number of filter stages: %d\n', getNumStages(D_multi_stage) ); fprintf('Stage 1: Decimation factor = %d , filter length = %d\n',... D_multi_stage.Stage1.DecimationFactor,... numel(D_multi_stage.Stage1.Numerator)); fprintf('Stage 2: Decimation factor = %d , filter length = %d\n',... D_multi_stage.Stage2.DecimationFactor,... numel(D_multi_stage.Stage2.Numerator)); fprintf('Stage 3: Decimation factor = %d , filter length = %d\n',... D_multi_stage.Stage3.DecimationFactor,... numel(D_multi_stage.Stage3.Numerator));
Number of filter stages: 3 Stage 1: Decimation factor = 6 , filter length = 33 Stage 2: Decimation factor = 2 , filter length = 11 Stage 3: Decimation factor = 4 , filter length = 101
Note that D_multi_stage is a three-stage multirate lowpass filter. We transform it to a bandpass filter by performing a frequency shift on the coefficients of each stage, while taking the cumulative decimation factor into account:
num1 = D_multi_stage.Stage1.Numerator; num2 = D_multi_stage.Stage2.Numerator; num3 = D_multi_stage.Stage3.Numerator; N1 = length(num1)-1; N2 = length(num2)-1; N3 = length(num3)-1; % Decimation factor between stage 1 & 2: M12 = D_multi_stage.Stage1.DecimationFactor; % Decimation factor between stage 2 & 3: M23 = D_multi_stage.Stage2.DecimationFactor; num1 = num1 .*exp(1j*(Fc / Fs)*2*pi*(0:N1)); num2 = num2 .*exp(1j*(Fc * M12 / Fs)*2*pi*(0:N2)); num3 = num3 .*exp(1j*(Fc * M12 * M23 / Fs)*2*pi*(0:N3)); D_multi_stage.Stage1.Numerator = num1; D_multi_stage.Stage2.Numerator = num2; D_multi_stage.Stage3.Numerator = num3;
Let's compare the cost of the single-stage and multistage filters:
fprintf('Single-stage filter cost:\n\n') cost(D_single_stage) fprintf('Multistage filter cost:\n') cost(D_multi_stage)
Single-stage filter cost: ans = struct with fields: NumCoefficients: 1129 NumStates: 1104 MultiplicationsPerInputSample: 23.5208 AdditionsPerInputSample: 23.5000 Multistage filter cost: ans = struct with fields: NumCoefficients: 113 NumStates: 140 MultiplicationsPerInputSample: 7.0208 AdditionsPerInputSample: 6.7500
Let's also compare the frequency response of the two filters:
vis = fvtool(D_single_stage,D_multi_stage,'DesignMask','off','legend','on'); legend(vis,'Single-stage','Multistage')
Let's use the multistage filter to perform zoom FFT:
fftlen = 32; L = BWFactor * fftlen; tones = [1625 2000 2125]; % sine wave tones sn = dsp.SineWave('SampleRate',Fs,'Frequency',tones,... 'SamplesPerFrame',L); Fsd = Fs / BWFactor; % Frequency points at which FFT is computed F = Fc + Fsd/fftlen*(0:fftlen-1)-Fsd/2; ap = dsp.ArrayPlot('XDataMode','Custom',... 'CustomXData',F,... 'YLabel','Magnitude',... 'XLabel','Frequency (Hz)',... 'YLimits',[0 1],... 'Title',sprintf('Zoom FFT. Resolution = %f Hz',(Fs/BWFactor)/fftlen)); for k=1:100 x = sum(sn(),2) + 1e-2 * randn(L,1); y = D_multi_stage(x); z = fft(y,fftlen,1); z = fftshift(z); ap( abs(z)/fftlen ) end release(ap)
dsp.ZoomFFT is a System object that implements zoom FFT based on the multirate multistage bandpass filter highlighted in the previous section. You specify the desired center frequency and decimation factor, and dsp.ZoomFFT will design the filter and apply it to the input signal.
Let's use dsp.ZoomFFT to zoom into the sine tones from the previous section's example:
zfft = dsp.ZoomFFT(BWFactor,Fc,Fs); for k=1:100 x = sum(sn(),2) + 1e-2 * randn(L,1); z = zfft(x); z = fftshift(z); ap( abs(z)/fftlen) end release(ap)
The zoom FFT block brings the functionality of dsp.ZoomFFT to Simulink. In the model dspzoomfft, we use the zoom FFT block to inspect the frequency band [800 Hz, 1600 Hz] of an input signal sampled at 44100 Hz.
 Multirate Signal Processing - Harris (Prentice Hall).
 Computing Translated Frequencies in digitizing and Downsampling Analog Bandpass - Lyons (https://www.dsprelated.com/showarticle/523.php)