FFT Analysis of Sine Sweep

42 Ansichten (letzte 30 Tage)
David Kendal
David Kendal am 23 Mai 2022
Kommentiert: Star Strider am 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');

Akzeptierte Antwort

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

Melden Sie sich an, um zu kommentieren.

Weitere Antworten (0)

Community Treasure Hunt

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

Start Hunting!

Translated by