MATLAB FFT bin size

I am attemping to create an FFT function that allows me to set the frequency range of the FFT, the number of bins, and windowing type. For the frequnecy range, my logic is to use a low-pass filter. If I add a filter, does it make sense to have windowing? How am I able to set a distinct number of bins in the FFT, it seems that automaticcaly the number of bins is incredibly high (I am collecting data at 48 kHz for about 90 seconds). I would like the bin width to be around 1 Hz instead of ~.01 which it currently is. How can I change this? We have another program that is able to achieve this, so I know that it is possible, but is this program cutting out data to achieve this?

[Steven Lord]

This sounds more like a general signal processing question than a question about MATLAB. If it is you might want to find a forum that is more targeted towards signal processing.

Or if you know what you want to do but aren't sure how to implement it in MATLAB, please show your code and ask a specific question about the section of the code where you're having doubts or uncertainty.

[dpb]

@Steven Lord is correct in that the subject is not specific to MATLAB at all, but I'll make one comment on a specific Q? raised...

"...I would like the bin width to be around 1 Hz instead of ~.01 which it currently is. How can I change this?"

The basic relationships are

Fmax = 1/(2*dt) Nyquist criterion

T = N*dt Sample data length

df = Fmax/(N/2) --> df = 2*Fmax/N Frequency bin width

Consequently, whatever you pick for N controls the length of the sampled data stream given a fixed sampling rate (dt). Your ability to discern frequency components in the signal without aliasing is totally dependent upon sampling at at least 2X the highest frequency component in the signal, but for reliable estimates, it really ought to be more like 3-4X. The only way to ensure a signal is not aliased is that the input either simply physically cannot contain any higher frequency data or by analog filtering prior to sampling--once the data are sampled, filtering has no effect; the aliasing has already occurred.

The FFT can be computed over any number of frequency bins up to the length of the input signal; select a number for it which will give you the desired frequency resolution. This will still process the whole time series but will aggregate the energy within each frequency bin; the total energy in the signal is always the same so half the number of bins will double the energy in each bin in order to integrate to the same total.

Answer 1

[Star Strider]

I agree with @Steven Lord, however there are ways to do what you want with the MATLAB fft function. I wrote a small utility function to do this (since a fair fraction of questions I choose to respond to involve calculating Fourier transforms). [SL: just fixed a few typos]

function [FTs1,Fv] = FFT1(s,t)

t = t(:);

L = numel(t);

if size(s,2) == L

s = s.';

end

Fs = 1/mean(diff(t));

Fn = Fs/2;

NFFT = 2^nextpow2(L);

FTs = fft((s - mean(s)) .* hann(L).*ones(1,size(s,2)), NFFT)/sum(hann(L));

Fv = Fs*(0:(NFFT/2))/NFFT;

% Fv = linspace(0, 1, NFFT/2+1)*Fn;

Iv = 1:numel(Fv);

FTs1 = FTs(Iv,:);

end

You can set the number of frequencies tto use with the second ‘NFFT’ argument. The frequencies themselves will go from 0 Hz (D-C) to the Nyquist frequency, however you can choose any frequencies less than those. Another option with respect to the frequencies is to use the interp1 function. I have never done this before, however using the ‘raw’ Fourier ttransform result my function returns, that would go something like this:

Freqs = linspace(0, Fn, 5E+4);

FTs1Intrp = interp1(Fv, FTs1, Feqs);

Testing that —

Fs = 4.8E+7;
Fn = Fs/2;
L = 1;                                                                  % One Second
t = linspace(0, Fs*L, Fs*L+1).'/Fs;                                     % Time (Column Vectors)
s = sum(sin(2*pi*(0:2.5E6:2.39E+7).*t), 2);                             % Signal Vector

[FTs1,Fv] = FFT1(s,t);

NrFreqs = numel(Fv)

NrFreqs = 33554433

figure
plot(Fv, abs(FTs1)*2)
grid
xlim('tight')
xlabel('Frequency (Hz)')
ylabel('Magnitude')
title('Original Frequencies')

Freqs = linspace(0, Fn, 5E+4);
FTs1Intrp = interp1(Fv, FTs1, Freqs);

figure
plot(Freqs, abs(FTs1Intrp)*2)
grid
xlim('tight')
xlabel('Frequency (Hz)')
ylabel('Magnitude')
title('Selected & Interpolated Frequencies')

function [FTs1,Fv] = FFT1(s,t)
    t = t(:);
    L = numel(t);
    if size(s,2) == L
        s = s.';
    end
    Fs = 1/mean(diff(t));
    Fn = Fs/2;
    NFFT = 2^nextpow2(L);
    FTs = fft((s - mean(s)) .* hann(L).*ones(1,size(s,2)), NFFT)/sum(hann(L));
    Fv = Fs*(0:(NFFT/2))/NFFT;
%     Fv = linspace(0, 1, NFFT/2+1)*Fn;
    Iv = 1:numel(Fv);
    FTs1 = FTs(Iv,:);
end

To decrease the maximum frequency displayed, no filtering is necessary. Just use the xlim function (with different arguments than I use here).

Also, 90 seconds at a 48 kHz sampling frequency is going to be a huge vector (4.32E6 elements or 34.56 MB).

.

[dpb]

"...90 seconds at a 48 kHz sampling frequency is going to be a huge vector."

And also raises the question of whether the signal is stationary or a transient -- doing baseband analysis on the whole thing will be smooshing whatever transient there is altogether.

As @Steven Lord hints, there's far more behind the question than a couple of quick responses about how MATLAB functions can be used can address...

[Star Strider]

Noted.

My intent is simply to demonstrate how to do the requested operations in MATLAB.

[dpb]

Noted. :)

My intent is to simply caution the OP that just because you can do something with MATLAB doesn't mean you'll get the correct result if the application of "that something" is not proper...

MATLAB is quite robust numerically but it can't tell you whether the assumptions the algorithms are based on are met or not. (Not that that is any news to you...)