5 views (last 30 days)

Show older comments

close all

clear all

clc

sine = [];

rate = 100000;

fbase=1300;

freq = fbase+100;

sinesine1 = 0.25*sin(linspace(0, 2*pi,rate/freq));

freq = fbase+200;

sinesine2 = 0.25*sin(linspace(0,2*pi,rate/freq));

freq = fbase+300;

sinesine3 =0.25*sin(linspace(0, 2*pi,rate/freq));

freq = fbase+400;

sinesine4 = 0.25*sin(linspace(0, 2*pi,rate/freq));

freq =fbase+500;

sinesine5 = 0.25*sin(linspace(0, 2*pi,rate/freq));

n=200;

for i = 1: n

sine = [sine sinesine1(1:end-1)];

end

for i = 1: n

sine = [sine sinesine2(1:end-1)];

end

for i = 1: n

sine = [sine sinesine3(1:end-1)];

end

for i = 1: n

sine = [sine sinesine4(1:end-1)];

end

for i = 1: n

sine = [sine sinesine5(1:end-1)];

end

Fs = rate;

L = length(sine); % Length of signal

T = 1/Fs; % Sampling period

t = (0:L-1)*T; % Time vector

Z = fft(sine);

P2= abs(Z)./L;

P1 = P2(1:L/2+1);

P1(2:end-1) = 2*P1(2:end-1);

f = Fs*(0:(L/2))/L;

figure('Name','Spettro 1400-1800','NumberTitle','off');

plot(f,P1)

axis([fbase fbase+1100 0 max(P5(2:end))+0.01*max(P5(2:end))])

Mathieu NOE
on 27 Apr 2021

hi again

I was not happy to see that your code needed a very high sampling rate to display correctly the 5 frequencies

this is by no mean required by Shannon's theorem

I modified the way the data are generated and now we can have the right display even with Fs = 5000;

also one single for loop needed vs five

close all

clearvars

clc

sine = [];

% rate = 100000;

% fbase= 1300;

% freq = fbase+100;

% sinesine1 = 0.25*sin(linspace(0, 2*pi,rate/freq));

% freq = fbase+200;

% sinesine2 = 0.25*sin(linspace(0,2*pi,rate/freq));

% freq = fbase+300;

% sinesine3 =0.25*sin(linspace(0, 2*pi,rate/freq));

% freq = fbase+400;

% sinesine4 = 0.25*sin(linspace(0, 2*pi,rate/freq));

% freq =fbase+500;

% sinesine5 = 0.25*sin(linspace(0, 2*pi,rate/freq));

% n=200;

% for i = 1: n

% sine = [sine sinesine1(1:end-1)];

% end

% for i = 1: n

% sine = [sine sinesine2(1:end-1)];

% end

% for i = 1: n

% sine = [sine sinesine3(1:end-1)];

% end

% for i = 1: n

% sine = [sine sinesine4(1:end-1)];

% end

% for i = 1: n

% sine = [sine sinesine5(1:end-1)];

% end

% Fs = rate;

% my suggestion

periods = 200;

Fs = 5e3;

dt = 1/Fs;

fbase= 1300;

freq = fbase + [100 200 300 400 500];

for ci = 1:length(freq)

omega = 2*pi*freq(ci);

t = (0:dt:periods/freq(ci)-dt);

sine = [sine 0.25*sin(omega*t)];

end

% T = 1/Fs; % Sampling period

% t = (0:L-1)*T; % Time vector

% L = length(sine); % Length of signal

% Z = fft(sine);

% P2= abs(Z)./L;

% P1 = P2(1:L/2+1);

% P1(2:end-1) = 2*P1(2:end-1);

% f = Fs*(0:(L/2))/L;

% figure('Name','Spettro 1400-1800','NumberTitle','off');

% plot(f,P1)

% axis([fbase fbase+1100 0 max(P5(2:end))+0.01*max(P5(2:end))])

%%%%%%%%%% PEAK HOLD AVERAGING FFT SPECTRUM DEMO %%%%%%%%%%%%%%%%%

NFFT = 250;

Overlap = 0.75;

[freq,fft_spectrum] = myfft_peak(sine(:), Fs, NFFT, Overlap); % data must be column oriented (samples x channels)

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

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

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

%%%%%%%%%% SPECTROGRAM DEMO %%%%%%%%%%%%%%%%%

[S,F,T] = myspecgram(sine, Fs, NFFT, Overlap);

figure(2);

imagesc(T,F,(abs(S)))

colorbar('vert');

set(gca,'YDir','Normal')

xlabel('Time (secs)')

ylabel('Freq (Hz)')

title('Short-time Fourier Transform spectrum')

colormap('jet');

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

% FFT peak spectrum of signal (example sinus amplitude 1 = 0 dB after fft).

% Linear averaging

% signal - input signal,

% Fs - Sampling frequency (Hz).

% nfft - FFT window size

% Overlap - buffer percentage of overlap % (between 0 and 0.95)

[samples,channels] = size(signal);

if channels > samples

signal = signal'; % flip dimensions

end

% fill signal with zeros if its length is lower than nfft

if samples<nfft

s_tmp = zeros(nfft,channels);

s_tmp((1:samples),:) = signal;

signal = s_tmp;

samples = nfft;

end

% window : hanning

window = hanning(nfft);

window = window(:);

% compute fft with overlap

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

spectnum = 1+ fix((samples-nfft)/offset); % Number of windows

% % for info is equivalent to :

% noverlap = Overlap*nfft;

% spectnum = fix((samples-noverlap)/(nfft-noverlap)); % Number of windows

% main loop

fft_spectrum = 0;

for i=1:spectnum

start = (i-1)*offset;

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

fft_spectrum = max(fft_spectrum,(abs(fft(sw))*4/nfft)); % X=fft(x.*hanning(N))*4/N; % hanning only

end

% fft_spectrum = fft_spectrum/spectnum; % to do linear averaging scaling

% one sidded fft spectrum % Select first half

if rem(nfft,2) % nfft odd

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

else

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

end

fft_spectrum = fft_spectrum(select,:);

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

end

function [fft_specgram,freq_vector,time] = myspecgram(signal, Fs, nfft, Overlap)

% FFT peak spectrogram of signal (example sinus amplitude 1 = 0 dB after fft).

% signal - input signal,

% Fs - Sampling frequency (Hz).

% nfft - FFT window size

% Overlap - buffer overlap % (between 0 and 0.95)

signal = signal(:);

samples = length(signal);

% fill signal with zeros if its length is lower than nfft

if samples<nfft

s_tmp = zeros(nfft,1);

s_tmp((1:samples),:) = signal;

signal = s_tmp;

samples = nfft;

end

% window : hanning

window = hanning(nfft);

window = window(:);

% compute fft with overlap

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

spectnum = 1+ fix((samples-nfft)/offset); % Number of windows

% % for info is equivalent to :

% noverlap = Overlap*nfft;

% spectnum = fix((samples-noverlap)/(nfft-noverlap)); % Number of windows

% main loop

fft_specgram = [];

for ci=1:spectnum

start = (ci-1)*offset;

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

fft_specgram = [fft_specgram abs(fft(sw))*4/nfft]; % X=fft(x.*hanning(N))*4/N; % hanning only

end

% one sidded fft spectrum % Select first half

if rem(nfft,2) % nfft odd

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

else

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

end

fft_specgram = fft_specgram(select,:);

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

% time vector

% time stamps are defined in the middle of the buffer

time = ((0:spectnum-1)*offset + round(offset/2))/Fs;

end

Mathieu NOE
on 26 Apr 2021

hello

FYI see code below

it iplement a peak hold fft averaging technique; also the spectrogram principle is coded (so you can see what matlab is doing)

both methods show the right amplitude

I reduced your sampling freq that was too high for the demo purpose

close all

clearvars

clc

sine = [];

rate = 10000;

fbase= 1300;

freq = fbase+100;

sinesine1 = 0.25*sin(linspace(0, 2*pi,rate/freq));

freq = fbase+200;

sinesine2 = 0.25*sin(linspace(0,2*pi,rate/freq));

freq = fbase+300;

sinesine3 =0.25*sin(linspace(0, 2*pi,rate/freq));

freq = fbase+400;

sinesine4 = 0.25*sin(linspace(0, 2*pi,rate/freq));

freq =fbase+500;

sinesine5 = 0.25*sin(linspace(0, 2*pi,rate/freq));

n=200;

for i = 1: n

sine = [sine sinesine1(1:end-1)];

end

for i = 1: n

sine = [sine sinesine2(1:end-1)];

end

for i = 1: n

sine = [sine sinesine3(1:end-1)];

end

for i = 1: n

sine = [sine sinesine4(1:end-1)];

end

for i = 1: n

sine = [sine sinesine5(1:end-1)];

end

Fs = rate;

% T = 1/Fs; % Sampling period

% t = (0:L-1)*T; % Time vector

% L = length(sine); % Length of signal

% Z = fft(sine);

% P2= abs(Z)./L;

% P1 = P2(1:L/2+1);

% P1(2:end-1) = 2*P1(2:end-1);

% f = Fs*(0:(L/2))/L;

% figure('Name','Spettro 1400-1800','NumberTitle','off');

% plot(f,P1)

% axis([fbase fbase+1100 0 max(P5(2:end))+0.01*max(P5(2:end))])

%%%%%%%%%% PEAK HOLD AVERAGING FFT SPECTRUM DEMO %%%%%%%%%%%%%%%%%

NFFT = 500;

Overlap = 0.75;

[freq,fft_spectrum] = myfft_peak(sine(:), Fs, NFFT, Overlap); % data must be column oriented (samples x channels)

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

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

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

%%%%%%%%%% SPECTROGRAM DEMO %%%%%%%%%%%%%%%%%

[S,F,T] = myspecgram(sine, Fs, NFFT, Overlap);

figure(2);

imagesc(T,F,(abs(S)))

colorbar('vert');

set(gca,'YDir','Normal')

xlabel('Time (secs)')

ylabel('Freq (Hz)')

title('Short-time Fourier Transform spectrum')

colormap('jet');

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

% FFT peak spectrum of signal (example sinus amplitude 1 = 0 dB after fft).

% Linear averaging

% signal - input signal,

% Fs - Sampling frequency (Hz).

% nfft - FFT window size

% Overlap - buffer percentage of overlap % (between 0 and 0.95)

[samples,channels] = size(signal);

if channels > samples

signal = signal'; % flip dimensions

end

% fill signal with zeros if its length is lower than nfft

if samples<nfft

s_tmp = zeros(nfft,channels);

s_tmp((1:samples),:) = signal;

signal = s_tmp;

samples = nfft;

end

% window : hanning

window = hanning(nfft);

window = window(:);

% compute fft with overlap

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

spectnum = 1+ fix((samples-nfft)/offset); % Number of windows

% % for info is equivalent to :

% noverlap = Overlap*nfft;

% spectnum = fix((samples-noverlap)/(nfft-noverlap)); % Number of windows

% main loop

fft_spectrum = 0;

for i=1:spectnum

start = (i-1)*offset;

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

fft_spectrum = max(fft_spectrum,(abs(fft(sw))*4/nfft)); % X=fft(x.*hanning(N))*4/N; % hanning only

end

% fft_spectrum = fft_spectrum/spectnum; % to do linear averaging scaling

% one sidded fft spectrum % Select first half

if rem(nfft,2) % nfft odd

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

else

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

end

fft_spectrum = fft_spectrum(select,:);

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

end

function [fft_specgram,freq_vector,time] = myspecgram(signal, Fs, nfft, Overlap)

% FFT peak spectrogram of signal (example sinus amplitude 1 = 0 dB after fft).

% signal - input signal,

% Fs - Sampling frequency (Hz).

% nfft - FFT window size

% Overlap - buffer overlap % (between 0 and 0.95)

signal = signal(:);

samples = length(signal);

% fill signal with zeros if its length is lower than nfft

if samples<nfft

s_tmp = zeros(nfft,1);

s_tmp((1:samples),:) = signal;

signal = s_tmp;

samples = nfft;

end

% window : hanning

window = hanning(nfft);

window = window(:);

% compute fft with overlap

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

spectnum = 1+ fix((samples-nfft)/offset); % Number of windows

% % for info is equivalent to :

% noverlap = Overlap*nfft;

% spectnum = fix((samples-noverlap)/(nfft-noverlap)); % Number of windows

% main loop

fft_specgram = [];

for ci=1:spectnum

start = (ci-1)*offset;

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

fft_specgram = [fft_specgram abs(fft(sw))*4/nfft]; % X=fft(x.*hanning(N))*4/N; % hanning only

end

% one sidded fft spectrum % Select first half

if rem(nfft,2) % nfft odd

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

else

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

end

fft_specgram = fft_specgram(select,:);

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

% time vector

% time stamps are defined in the middle of the buffer

time = ((0:spectnum-1)*offset + round(offset/2))/Fs;

end

Mathieu NOE
on 27 Apr 2021

Ok

so let's put back the rate at 100 kHz, then the 5 peaks will show up but you'll have to increase the NFFT to 10 times higher if you'd like to distinguish the relatively close frequencies

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!