Decreasing bin size FFT options for increased resolution

I have a data set sampled at 0.8Mhz that is an impact blast wave that has 6ms capture leading to about 5000 points of data. We are interested in doing an FFT to look at the frequency of the data set. However, due to the sameple rate and sameple length we only end up with 2500 points meaning our bins are ~160Hz bins, which does not carry enough data for our analysis we are interested in. Our team was exploring ways to increase the data content of the wave so we can increase the granulariuty of the frequency domain.

Our first idea which we wanted to check was if it was "ok" to copy and paste the data along the X axis in time so that we would not add any spectral leckage into the system while decrasing the bin size for increased resolution. I did not see many comments on this so I am not sure if this is taboo or not.

There is the option of zero padding but that leads to spectral leakage and due to the nature of this data, we are unsure if this is viable with the sinc(x) function being used for this.

The end goal is a PSD using pwelch or periodogram with a force-exponential window followed by a bode plot of magnitude and phase ( which I am also only seeing matlab functions to be used on transfer system functions and havnt found a way to perform on a raw data set.)

So I am looking for insights on how best to aproach this.

A follow up question that I was having a hard time understanding in my signals and systems textbooks that i have been harvesting info out of is whenyou do an fft, the data point that has the stem point ( in this case the first woud be 160Hz) is that the data in that exact frequency or is that a combination of all of the frequencies not reported after performing the transformation with bleed between the points.

Thank you for the help!

Answer 1

[Star Strider]

You can minimise the spectral leakage by windowing the fft (except for a rectangular window, which is the default if no others are chosen). Subtracting the mean from the signal before calculating the fft elimniates the D-C offset, making the other peaks more easily visible.

Most of the information is below 5 (Hz?), so I limited the frequency-domain plot to that range.

Try this —

LD = load('datainquestion.mat');
a = LD.a;
t = a(:,1) * 1E+3;                                              % Time (Converted To milliseconds)
psi = a(:,2);

figure
plot(t, psi)
grid
xlabel('Time (ms)')
ylabel('Amplitude (units)')

L = size(a,1);                                                  % Signal Length
Ts = mean(diff(t));                                             % Sampling Interval
% Tsd = std(diff(x))                                              % Check Sampling Interval Variation
Fs = 1/Ts                                                       % Sampling Frequency

Fs = 0.8000

Fn = Fs/2;                                                      % Nyquist Frequency

NFFT = 2^(nextpow2(L)+6)                                        % Length Of 'fft' (Integer Power-Of-2 Is Most Efficient)

NFFT = 524288

hw = hann(L);                                                   % Window Function
FTpsi = fft((psi - mean(psi)).*hw, NFFT)/sum(hw);               % Calculate Windowed 'fft', Subtract Mean To Eliminate D-C Offset
Fv = Fs*(0:(NFFT/2))/NFFT;                                      % Frequency Vector
Iv = 1:numel(Fv);                                               % Index Vector

figure
plot(Fv, abs(FTpsi(Iv))*2)
grid
xlim([0 2.55E-3])
xlabel('Frequency (Hz)')
ylabel('Magnitude')

figure
plot(Fv, mag2db(abs(FTpsi(Iv))*2))
grid
xlabel('Frequency (Hz)')
ylabel('Power (dB)')
xlim([0 2.55E-3])

I’m not certain that this meets your requirements. It’s how I usually calculate the fft. (If you want the maximum to be at 0 dB, divide ‘abs(FTpsi))’ by ‘max(abs(FTpsi))’.)

EDIT — (19 Apr 2024 at 18:40)

Converted ‘x’ to ‘t’ and converted ‘t’ to milliseconds (multiplying it by 1E+3).

EDIT — (19 Apr 2024 at 19:00)

Changed the xlim range to [0 2.55E-3].

.

[Kyle]

Sweet. You got to this really fast!

But I do have some questions. 1st the time I posted, should have been labeled in ms. Because of that the Fs was only at 800 for the calcuations and not 800,000. Because we are working in a different time scale, and the frequency is dependent on the sampling rate, I assume it is important detail that I overlooked in my explaniation as the NFFT is only 16384.

In the way that you plotted it, is there any way to increase the quantaity of points between [0-2] which was the bulk of the question as you are spot on. The large % of the content is below 5KHz, we were interested in focusing our analysis on the 2KHz band but as we only have so few data points we found that challeneging.

2nd why did you use the mag2db function instead of a pwlech or a specogram function. I would tend to think all 3 do a similar task?

Thank you for the help as always.

[Star Strider]

Thank you!

1. I just now converted the time units to milliseconds by multiplying them by .

The value for ‘NFFT’ can be arbitrarily large, depending on what you want, however it is always best to have it as an integer power-of-2. Indreasing it to:

NFFT = 2^(nextpow2(L)+6)

(equalling 524288) does not significantly improve the frequency resolution (visually) beyond just adding 2 to it instead of 6. It definitely increases the frequency resolution.

2. I used mag2db to convert the magnitude to power using the simple fft result. The other functions you mention, pwelch and spectrogram, calculate and plot the power as power spectral density, or dB/Hz. My plot data are simply dB without converting to spectral density. (This is the same result that the pspectrum funciton would provide.) I find this more intuitive, so simply a personal preference. You are of course free to use whatever functions you want.

My pleasure!

.