1
$\begingroup$

I am trying to calculate the area under my FIR filter kernel using trapz but I'm getting weird results. The reason I want to know the area is eventually I would like to loop through filter orders in a small range and pick the order which most closely matches my ideal filter window.

When you run this code in Matlab the area under the kernel is given as 9.5e3. However, if you look at my ideal filter window its area is 2e3.

What am I missing here and why is my integral result way off?

Matlab Code:

clc,clear;
close all;

srate = 250e3;
nyquist = srate/2;
npnts = srate*5;    %generate number of points for 5 seconds of sampling
time = (0:npnts-1)/srate; 
hz = linspace(0,srate/2,floor(npnts/2)+1)/1e3;  %go from 0 to nyquist, frequency resolution is defined by last term


% FIR filter specs
passband = [25e3 27e3];
transw = 0.01;   %how the filter's edges taper
test_order = 13;

order = round(test_order*srate/passband(1)); % define the number of time points for the filter kernel

shape = [ 0 0 1 1 0 0 ];    %FIR shape

%define frequency shape of the FIR shape, firls function requires
%frequencies to by normalized by nyquist
frex = [0 passband(1)-passband(1)*transw passband passband(2)+passband(2)*transw nyquist]/nyquist;

filter_kernel = firls(order, frex, shape);

% power spectrum of filter kernel
filtkernX = abs(fft(filter_kernel,npnts)).^2;

% calculate area under filter kernel 
disp(trapz( filtkernX(1:length(hz))) )

%plot
figure;
plot(hz,filtkernX(1:length(hz)),'-','linew',2,'markersize',1)
hold on 
plot([0 passband(1) passband passband(2) nyquist]./1e3,[0 0 1 1 0 0],'ro-','linew',2,'markerfacecolor','w')
xlim([0 60])
xlabel('Frequency (kHz)')
ylabel('Amplitude a.u.')
title('Frequency Response')
legend({'Actual';'Ideal'})

$\endgroup$
0

1 Answer 1

3
$\begingroup$

However, if you look at my ideal filter window its area is 2e3.

Nope. It's 10000. The spectrum is a metric of spectral density, not frequency. In your case your bandwidth is 2000 Hz and the frequency resolution is 0.2Hz ($\Delta f = \frac{f_s}{N_{FFT}}$). That means your passband has $M = \frac{2000Hz}{0.2Hz} = 10000$ points in the DFT and if you integrate 10000 ones you get 10000.

Side comment: in most cases it's preferable to integrate the power, not the amplitude.

$\endgroup$
2
  • $\begingroup$ If I were to switch to integrating over power, would the integral of the ideal filter window stay the same? $\endgroup$ Apr 6, 2023 at 19:03
  • 1
    $\begingroup$ Yes. The square of 1 is 1. $\endgroup$
    – Hilmar
    Apr 6, 2023 at 19:17

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.