FFT Analysis of Sine Sweep

34 visualizaciones (últimos 30 días)
David Kendal
David Kendal el 23 de Mayo de 2022
Comentada: Star Strider el 24 de Mayo de 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');

Respuesta aceptada

Star Strider
Star Strider el 23 de Mayo de 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 comentarios
David Kendal
David Kendal el 24 de Mayo de 2022
Also for these I notice the units of power spectral density are both different 🤷🏻♂ cheers
Star Strider
Star Strider el 24 de Mayo de 2022
They may be calculated differently, however I doubt that hte actual values are much different.

Iniciar sesión para comentar.

Más respuestas (0)

Community Treasure Hunt

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

Start Hunting!

Translated by