clc

data = importdata('SFC5S_nov25_ST_1_1.txt');

dt = 1e-6;

signal = data(:,1);

samples = length(signal);

Fs = 1/dt;

decim = 1;

if decim>1

signal = decimate(signal,decim);

Fs = Fs/decim;

end

NFFT = 512;

OVERLAP = 0.95;

spectrogram_dB_scale = 80;

option_w = 0;

[freq, sensor_spectrum] = myfft_peak(signal,Fs,NFFT,OVERLAP);

sensor_spectrum_dB = 20*log10(sensor_spectrum);

if option_w == 1

pondA_dB = pondA_function(freq);

sensor_spectrum_dB = sensor_spectrum_dB+pondA_dB;

my_ylabel = ('Amplitude (dB (A))');

else

my_ylabel = ('Amplitude (dB (L))');

end

figure(1),plot(freq,sensor_spectrum_dB,'b');grid

title(['Averaged FFT Spectrum / Fs = ' num2str(Fs) ' Hz / Delta f = ' num2str(freq(2)-freq(1)) ' Hz ']);

xlabel('Frequency (Hz)');ylabel(my_ylabel);

[sg,fsg,tsg] = specgram(signal,NFFT,Fs,hanning(NFFT),floor(NFFT*OVERLAP));

sg_dBpeak = 20*log10(abs(sg))+20*log10(2/length(fsg));

if option_w == 1

pondA_dB = pondA_function(fsg);

sg_dBpeak = sg_dBpeak+(pondA_dB*ones(1,size(sg_dBpeak,2)));

my_title = ('Spectrogram (dB (A))');

else

my_title = ('Spectrogram (dB (L))');

end

min_disp_dB = round(max(max(sg_dBpeak))) - spectrogram_dB_scale;

sg_dBpeak(sg_dBpeak<min_disp_dB) = min_disp_dB;

figure(2);

imagesc(tsg,fsg,sg_dBpeak);colormap('jet');

axis('xy');colorbar('vert');grid

title([my_title ' / Fs = ' num2str(Fs) ' Hz / Delta f = ' num2str(fsg(2)-fsg(1)) ' Hz ']);

xlabel('Time (s)');ylabel('Frequency (Hz)');

function pondA_dB = pondA_function(f)

n = ((12200^2*f.^4)./((f.^2+20.6^2).*(f.^2+12200^2).*sqrt(f.^2+107.7^2).*sqrt(f.^2+737.9^2)));

r = ((12200^2*1000.^4)./((1000.^2+20.6^2).*(1000.^2+12200^2).*sqrt(1000.^2+107.7^2).*sqrt(1000.^2+737.9^2))) * ones(size(f));

pondA = n./r;

pondA_dB = 20*log10(pondA(:));

end

function [freq_vector,fft_spectrum] = myfft_peak(signal, Fs, nfft, Overlap)

signal = signal(:);

samples = length(signal);

if samples<nfft

s_tmp = zeros(nfft,1);

s_tmp((1:samples)) = signal;

signal = s_tmp;

samples = nfft;

end

window = hanning(nfft);

window = window(:);

offset = fix((1-Overlap)*nfft);

spectnum = 1+ fix((samples-nfft)/offset);

fft_spectrum = 0;

for i=1:spectnum

start = (i-1)*offset;

sw = signal((1+start):(start+nfft)).*window;

fft_spectrum = fft_spectrum + (abs(fft(sw))*4/nfft);

end

fft_spectrum = fft_spectrum/spectnum;

if rem(nfft,2)

select = (1:(nfft+1)/2)';

else

select = (1:nfft/2+1)';

end

fft_spectrum = fft_spectrum(select);

freq_vector = (select - 1)*Fs/nfft;

end

## Best Answer