FFT Analysis of Sine Sweep

62 vues (au cours des 30 derniers jours)
David Kendal
David Kendal le 23 Mai 2022
Commenté : Star Strider le 24 Mai 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');

Réponse acceptée

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

Connectez-vous pour commenter.

Plus de réponses (0)

Community Treasure Hunt

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

Start Hunting!

Translated by