FFT Analysis of Sine Sweep

42 views (last 30 days)
David Kendal
David Kendal on 23 May 2022
Commented: Star Strider on 24 May 2022
Hi there, I am trying to run FFT analysis of the sine sweep I have created and not sure if the graphs I'm producing are working as they should. I feel they look to smooth and wanted to see if I could do anything better. The code runs and works but I still feel there is something missing. Any help greatly appreciated.
Many thanks
T=5; %size of window
fs=44100; %sampling frequency
df=1/T; %frequency res
dt=1/fs; %time resolution
t=(0:+dt:T-dt); %time vector
df_t=500; %swept rate (Hz/seconds)
% pre-allocate size of t:
sweptsin = zeros(size(t));
for i=1:+1:length(t)
%i=i+1;
if(i==1) %initialise f and t.
f=20; ti=0;
else
ti=ti+dt; %time increment
f=f+df_t*dt; %freq increment
end
w=2*pi*f; %omega
sweptsin(i)=sin(w*ti); %swept sine wave
end
NFFT=1024; %NFFT-point DFT
X=fftshift(fft(sweptsin,NFFT)); %compute DFT using FFT
fVals1=(-NFFT/2:NFFT/2-1)/NFFT; %DFT Sample points
subplot(4,1,1)
plot(fVals1,abs(X));
title('Double Sided FFT - with FFTShift');
xlabel('Normalized Frequency')
ylabel('DFT Values');
L=length(sweptsin);
X=fft(sweptsin,NFFT);
Px=X.*conj(X)/(NFFT*L); %Power of each freq components
fVals2=fs*(0:NFFT/2-1)/NFFT;
subplot(4,1,2)
plot(fVals2,Px(1:NFFT/2),'b','LineSmoothing','on','LineWidth',1);
title('One Sided Power Spectral Density');
xlabel('Frequency (Hz)')
ylabel('PSD');
X=fft(sweptsin,NFFT); %compute DFT using FFT
nVals=(0:NFFT-1)/NFFT; %Normalized DFT Sample points
subplot(4,1,3)
plot(nVals,abs(X));
title('Double Sided FFT - without FFTShift');
xlabel('Normalized Frequency')
ylabel('DFT Values');
NFFT=1024;
L=length(sweptsin);
X=fftshift(fft(sweptsin,NFFT));
Px=X.*conj(X)/(NFFT*L); %Power of each freq components
fVals3=fs*(-NFFT/2:NFFT/2-1)/NFFT;
subplot(4,1,4)
plot(fVals3,10*log10(Px),'b');
title('Power Spectral Density');
xlabel('Frequency (Hz)')
ylabel('Power');

Accepted Answer

Star Strider
Star Strider on 23 May 2022
Your posted code is likely as good as it is possible to get.
I added a pspectrum call at the end —
T=5; %size of window
fs=44100; %sampling frequency
df=1/T; %frequency res
dt=1/fs; %time resolution
t=(0:+dt:T-dt); %time vector
df_t=500; %swept rate (Hz/seconds)
% pre-allocate size of t:
sweptsin = zeros(size(t));
for i=1:+1:length(t)
%i=i+1;
if(i==1) %initialise f and t.
f=20; ti=0;
else
ti=ti+dt; %time increment
f=f+df_t*dt; %freq increment
end
w=2*pi*f; %omega
sweptsin(i)=sin(w*ti); %swept sine wave
end
NFFT=1024; %NFFT-point DFT
X=fftshift(fft(sweptsin,NFFT)); %compute DFT using FFT
fVals1=(-NFFT/2:NFFT/2-1)/NFFT; %DFT Sample points
subplot(4,1,1)
plot(fVals1,abs(X));
title('Double Sided FFT - with FFTShift');
xlabel('Normalized Frequency')
ylabel('DFT Values');
L=length(sweptsin);
X=fft(sweptsin,NFFT);
Px=X.*conj(X)/(NFFT*L); %Power of each freq components
fVals2=fs*(0:NFFT/2-1)/NFFT;
subplot(4,1,2)
plot(fVals2,Px(1:NFFT/2),'b','LineSmoothing','on','LineWidth',1);
Warning: The LineSmoothing property will be removed in a future release.
Warning: The LineSmoothing property will be removed in a future release.
title('One Sided Power Spectral Density');
xlabel('Frequency (Hz)')
ylabel('PSD');
X=fft(sweptsin,NFFT); %compute DFT using FFT
nVals=(0:NFFT-1)/NFFT; %Normalized DFT Sample points
subplot(4,1,3)
plot(nVals,abs(X));
title('Double Sided FFT - without FFTShift');
xlabel('Normalized Frequency')
ylabel('DFT Values');
NFFT=1024;
L=length(sweptsin);
X=fftshift(fft(sweptsin,NFFT));
Px=X.*conj(X)/(NFFT*L); %Power of each freq components
fVals3=fs*(-NFFT/2:NFFT/2-1)/NFFT;
subplot(4,1,4)
plot(fVals3,10*log10(Px),'b');
title('Power Spectral Density');
xlabel('Frequency (Hz)')
ylabel('Power');
figure
pspectrum(sweptsin,t,'spectrogram')
colormap(turbo)
.
  4 Comments
David Kendal
David Kendal on 24 May 2022
Also for these I notice the units of power spectral density are both different 🤷🏻♂ cheers
Star Strider
Star Strider on 24 May 2022
They may be calculated differently, however I doubt that hte actual values are much different.

Sign in to comment.

More Answers (0)

Community Treasure Hunt

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

Start Hunting!