Spectogram ve Chirp Sinyali Açıklamasıyla Matlab Kodları

Spectogram Nedir? Spectogram aslında sinyalin parça parça DTFT alarak frekansın zamanla değişimi gözlemlemeye yarar.

Sinyali parçalara ayırmak için zaman düzleminde pencereleme işlemi yapıyor. Matlab Spectogram varsayılan pencereleme fonksiyonu Hamming algoritmasıdır. * s = spectrogram(x,window,noverlap,nfft); * x = girdi sinyal * window= blackman(128) yada hamming(256) pencereleme algoritması ve zamanı parçalama sayısı. Spectogram çözünürlüğü bu sayı arttıkça artar. 2 ve üssü seçmekte yarar var. matlab spectogram içinde fft kullanıyor ve fft 2 ve üssü ile hesap yapar. Sinyal (1-cos(2*pi*(0:N)'/N))/2 N parçadan oluşuyorsa hann penceresi N+1 uzunluğundadır. * noverlap= positif tam sayıdır. Örneklerin üst üste gelme sayısı. Parçalar arasında varsayılan olarak %50 olacak şekilde üst üste biner. * nfft= frekans düzleminin örnekleme sayısı. Genelde girdi sinyalinin örnekleme frekansına eşit olur.

Contents

• Most Known Windowing Techniques:
• Kaiser, or Kaiser-Bessel, window:
• Rectangular(hamming) window: this is really no window, i.e. no tapering. It is mentioned for completeness and because the rectangular window is one of the window options in many window-related VIs in Labview.
• Hangisini Seçmeliyiz?
• Matlab Window Designer
• Yazdığım spectogram. Matlab kadar iyi değil ama stftnin fft ile yapıldığını kanıtlar.
• Spectogram can invertable. Inverse STFT
• FMCW
• Matlab chirp sinyali.
• Logaritmik chirp
• Kendi chirp sinyalin...
• spectrum
• Temel Linear FMCW Modulasyonu
• 1.Linear FMCW modulasyon testi.
• Matlab Radar waveform analyzer ürettiği sinyal
• Close and clear everything

Most Known Windowing Techniques:

1. Hanning window: This is also known as a cosine taper. It starts at 0, rises to 1 in the middle of the period, and then goes smoothly back down to zero at the end. w1(t) = 0.5 – 0.5 * cos(2pt/T) If x(t) is the data (presumed to have zero mean value), then the windowed version of the data,y(t), is computed as follows: y1(t) = x(t) * w1(t)

1. Hamming window: This modified cosine taper starts at 0.08, rises to 1 in the middle of the period, and then goes smoothly back down to 0.08 at the end. w2(t) = 0.54 – 0.46 * cos(2pt/T) y2(t) = x(t) * w2(t)

1. Blackman window: This modified cosine taper uses two cosines. w3(t) * 0.42 – 0.5 * cos(2pt/T) + 0.08 * cos(4pt/T) y3(t) = x(t) * w3(t)

w1=hanning(100);
w2=hamming(100);
w3=blackman(100);
t=1:100;
plot(t,w1,'.r-',t,w2,'.g-',t,w3,'.b-');
xlabel('Time'); ylabel('Amplitude'); title('Window Comparison');
legend('Hanning','Hamming','Blackman'); grid on;

Kaiser, or Kaiser-Bessel, window:

This tapering function has an adjustable parameter b. b can vary from 0 to +infinity, but the range 3<=b<=10 is typical. When b = 0, the Kaiser window becomes a rectangular window (no tapering). As b rises, tapering increases, and the main lobe width increases and sideband attenuation increases. The exact formula for the Kaiser window is complicated. Figure 2 shows Kaiser windows with b=3, 5, and 10.

w4=kaiser(100,3);
w5=kaiser(100,5);
w6=kaiser(100,10);
t=1:100;
plot(t,w4,'.r-',t,w5,'.g-',t,w6,'.b-');
xlabel('Time'); ylabel('Amplitude'); title('Kaiser Window Comparison');
legend('b=3','b=5','b=10'); grid on;

Rectangular(hamming) window: this is really no window, i.e. no tapering. It is mentioned for completeness and because the rectangular window is one of the window options in many window-related VIs in Labview.

w0(t)=1 for 0<=t<T. y(t) = x(t)*w0(t) = x(t)

t=0:1/1e2:2;
d=chirp(t,0,1,250);
N=length(d);
x=d-mean(d); % sıfır frekans bileşeni silinir.
w1=hamming(N);
y=x.*w1;
plot(t,x,'k',t,y,'r');
xlabel('Time (s)'); ylabel('Amplitude (V)'); title('chirp');
legend('Raw (zero mean)','After Hamming Window');

Hangisini Seçmeliyiz?

Genelde yan lopların büyük olmaması, analobun olabildiğince dar olması, sisteme yeni harmonik katmaması ve filtre etmek istemediğiniz bölgelere en az etki yapması tercih sebebidir. Fakat bu özellikleri ideal olarak sağlamak imkansıza yakındır. analob daraldıkça işlem gücüde orantılı olarak artar.

w0=ones(1,100);
wvtool(w1,w2,w3,w0);

Matlab Window Designer

filt=load('filtreler.mat');
filtAr=struct2cell(filt).';
N=13;
tit=['hamming64','hann64','rectengular64','triangular64','barlett64','barletHannig64','blackman64','blackmanHarris64',...
    'chebyshev64','gauusian64','kaiser64','kaiser32','kaiser16'];
disp('pencelerin görünümleri ve frekans tepkileri.')
figure('Name','Time Domain')

%Time Domain
for p=1:N
    subplot(4,N/4+0.75,p)
    plot(filtAr{p})
    xlabel('samples')
    ylabel('amplitude')
    title(tit(p))
%     H=wvtool(filtAr{p});
%     H.Name=tit(p);
end

% Frequency Domain
figure('Name','Frequency Domain Response')
for p=1:N
    subplot(4,N/4+0.75,p)

    [H,W] = freqz(filtAr{p});
    plot(W,20*log10(abs(H)))

    xlabel('Normalized Frequency x pi rad/sample')
    ylabel('Magnitude (dB)')
    title(tit(p))
%     H=wvtool(filtAr{p});
%     H.Name=tit(p);
end

pencelerin görünümleri ve frekans tepkileri.

Yazdığım spectogram. Matlab kadar iyi değil ama stftnin fft ile yapıldığını kanıtlar.

Teori: Sinyali eşit belli aralıklara bölerek her birinin ayrı ayrı fft alınmasıdır En çok müzik analizi yada radar sinyal analizi yapılırken kullanılır.

frekans değişimi zaman değişimiyle varsayılan olarak y ekseninde görülebilir.

clear, clc, close all

fs=1e3;
tstart=0;
tstop=2;
t = tstart:1/fs:tstop;
x = chirp(t,0,1,250); % t=0 => f0=0 Hz, t=1 => f1=250 Hz

% define analysis parameters
xlen = length(x);                   % length of the signal
wlen = 1024;                        % window length (recomended to be power of 2)
hop = wlen/4;                       % hop size (recomended to be power of 2)
nfft = 4096;                        % number of fft points (recomended to be power of 2)

% perform STFT
[S, f, t] = stft(x, wlen, hop, nfft, fs);

% define the coherent amplification of the window
K = sum(hamming(wlen, 'periodic'))/wlen;

% take the amplitude of fft(x) and scale it, so not to be a
% function of the length of the window and its coherent amplification
S = abs(S)/wlen/K;

% correction of the DC & Nyquist component
if rem(nfft, 2)                     % odd nfft excludes Nyquist point
    S(2:end, :) = S(2:end, :).*2;
else                                % even nfft includes Nyquist point
    S(2:end-1, :) = S(2:end-1, :).*2;
end

% convert amplitude spectrum to dB (min = -120 dB)
S = 20*log10(S + 1e-6);

% plot the spectrogram
figure(1)
surf(t, f, S)
shading interp
axis tight
box on
view(0, 90)
set(gca, 'FontName', 'Times New Roman', 'FontSize', 14)
xlabel('Time, s')
ylabel('Frequency, Hz')
title('Amplitude spectrogram of the signal')

handl = colorbar;
set(handl, 'FontName', 'Times New Roman', 'FontSize', 14)
ylabel(handl, 'Magnitude, dB')

Spectogram can invertable. Inverse STFT

clear, clc, close all

fs=1e3;
tstart=0;
tstop=2;
t = tstart:1/fs:tstop;
x = chirp(t,0,1,250); % t=0 => f0=0 Hz, t=1 => f1=250 Hz

% signal parameters
xlen = length(x);
t = (0:xlen-1)/fs;

% define the analysis and synthesis parameters
wlen = 1024;
hop = wlen/4;
nfft = 10*wlen;

% perform time-frequency analysis and resynthesis of the original signal
[stft, f, t_stft] = stft(x, wlen, hop, nfft, fs);
[x_istft, t_istft] = istft(stft, wlen, hop, nfft, fs);

% plot the original signal
figure(1)
plot(t, x, 'b')
grid on
xlim([0 max(t)])
set(gca, 'FontName', 'Times New Roman', 'FontSize', 14)
xlabel('Time, s')
ylabel('Signal amplitude')
title('Original and reconstructed signal')

% plot the resynthesized signal
hold on
plot(t_istft, x_istft, '-.r')
legend('Original signal', 'Reconstructed signal')

FMCW

Projede sürekli dalga modulasyonu yapabilmeyi kolaylaştıracak bir fonksiyon yazacağım. **Bu modulasyon sayesinde frekans-zaman grafiğinde içine konan bir fonksiyonun görüneceği gibi frekans değişimi sağlanacaktır. *Eğer yapabilirsem aynı zamanda verilen fonksiyona göre genlik değişimini gösterebilirim. s(t)=A*cos(2*pi*f0*t) input=fonksiyon output=frekansı değişmiş sürekli dalga

Matlab chirp sinyali.

Linear chirp y=chirp(t,f0,t1,f1); y = chirp(t,f0,t1,f1,'quadratic',phi,'shape');

fs=1e3;
tstart=0;
tstop=2;
t = tstart:1/fs:tstop;
f0=0;
t1=tstop; f1=250; %Hz
y = chirp(t,f0,t1,f1); % t=0 => f0=0 Hz, t=1 => f1=250 Hz
figure;
subplot(2,1,1)
plot(t,y)
subplot(2,1,2)
spectrogram(y,256,250,256,fs,'yaxis') %sampling frequency must be same or high than signal sampling frequency

Logaritmik chirp

fs=1e3;
tstart=0;
tstop=2;
t = tstart:1/fs:tstop;
f0=10;
t1=tstop; f1=250; %Hz
y = chirp(t,f0,t1,f1,'logarithmic'); % t=0 => f0=0 Hz, t=1 => f1=250 Hz
figure;
subplot(2,1,1)
plot(t,y)
subplot(2,1,2)
spectrogram(y,256,250,256,fs,'yaxis') %sampling frequency must be same or high than signal sampling frequency

Kendi chirp sinyalin...

fs=500; %sampling frequency
t=0:1/fs:1; %time base - upto 1 second

f0=1;% starting frequency of the chirp
f1=fs/20; %frequency of the chirp at t1=1 second
x = mychirp(t,f0,1,f1);
subplot(2,2,1)
plot(t,x,'k');
title('Chirp Signal');
xlabel('Time(s)');
ylabel('Amplitude');

L=length(x);
NFFT = 1024;
X = fftshift(fft(x,NFFT));
Pxx=X.*conj(X)/(NFFT*NFFT); %computing power with proper scaling
f = fs*(-NFFT/2:NFFT/2-1)/NFFT; %Frequency Vector

subplot(2,2,2)
plot(f,abs(X)/(L),'r');
title('Magnitude of FFT');
xlabel('Frequency (Hz)')
ylabel('Magnitude |X(f)|');
xlim([-50 50])
Pxx=X.*conj(X)/(NFFT*NFFT); %computing power with proper scaling
subplot(2,2,3)
plot(f,10*log10(Pxx),'r');
title('Double Sided - Power Spectral Density');
xlabel('Frequency (Hz)')
ylabel('Power Spectral Density- P_{xx} dB/Hz');
xlim([-100 100])

X = fft(x,NFFT);
X = X(1:NFFT/2+1);%Throw the samples after NFFT/2 for single sided plot
Pxx=X.*conj(X)/(NFFT*NFFT);
f = fs*(0:NFFT/2)/NFFT; %Frequency Vector
subplot(2,2,4)
plot(f,10*log10(Pxx),'r');
title('Single Sided - Power Spectral Density');
xlabel('Frequency (Hz)')
ylabel('Power Spectral Density- P_{xx} dB/Hz');

spectrum

figure
[S,F,T,p] = spectrogram(x,blackman(256),[],128,fs);
[q,nd] = max(10*log10(p));
% % helperCWTTimeFreqPlot(S,T,F,'surf','STFT of Quadratic Chirp','Seconds','Hz')
surf(T, F, abs(S))
hold on
plot3(T,F(nd),q,'r','linewidth',4)
hold off

shading interp
axis tight
box on
view(0, 90) % view(2)
set(gca, 'FontName', 'Times New Roman', 'FontSize', 14)
xlabel('Time, s')
ylabel('Frequency, Hz')
zlabel('Amplitude V')
title('Amplitude spectrogram of the signal')

handl = colorbar;
set(handl, 'FontName', 'Times New Roman', 'FontSize', 14)
ylabel(handl, 'Magnitude, dB')

Temel Linear FMCW Modulasyonu

Faz fit içine eklenebilir ama faz değerini 0 olarak tuttum. *Burada amacım sinuzoidal fonksiyonun içine istediğiniz fonksiyonu yerleştirip istediğimiz gibi chirp fonksiyonu elde etmek. Burada beta yani K=beta=(f1?f0)/t1 olarak sabit tuttum. Fakat beta(K) katsayısı frekans değişim fonksiyonuyla hep aynı hesaplanmadığından istediğimiz sinyali tam elde edemiyoruz. Teori; fi(t)=f0+?*t Lineer frekans değişim fonksiyonumuzdur. Buradan ?=(f1?f0)/t1 şeklinde hesaplanır. fi(t)=f0+?*t.^2 Eğer t.^2 frekans değişimiyle yapmak istersek ?=(f1?f0)/(t1.^2) olarak hesaplanır. fi(t)=f0×?.^t spectrumda logaritmik bir ifade üretir. ?=(f1f0).^(1/t1) ; [f1 > f0] olarak hesaplanır. K pozitif ise frekansı arttırır. veya negatif ise frekansı azaltır.

modFreqFonk=@(x,offset) (x+offset);
fit=@(f0,t,K,x,offset) 2*pi*(f0+(K/2).*modFreqFonk(x,offset)).*t; % -T/2<= t <=T/2
st=@(A,f0,t,K,x,offset) A.*exp(1i*fit(f0,t,K,x,offset)); %real kısmı cos oluyor.

% Test Verileri
f0=10;
tstart=0;
tstop=2;
A=1;

t0=tstart;
t1=tstop;
T=t1-t0;
f11=f0;
f00=1;
K=(f11-f00)/T;

fskat=40;
fs=f0.*abs(K.*fskat); %frekans artışına karşın örnekleme frekansının artışı.
if K==0
    fs=f0.*abs(fskat);
end
% fs=1e3;

t=tstart:1/fs:tstop;
x=(linspace(-10,10,numel(t)) );
offset=0;

% L=length(t); % Sinyalin uzunluğu
divi=2;

1.Linear FMCW modulasyon testi.

zaman-frekans grafiği düzgün değil.

tic;
mychirp2(modFreqFonk,st,A,x,offset,f0,fs,fskat,t,K,divi); toc

Elapsed time is 0.122183 seconds.

Matlab Radar waveform analyzer ürettiği sinyal

chirpTriangle

Close and clear everything

close all;clear;clc;

Published with MATLAB® R2017b

Bu dökümanı "matlab publish" ile ürettim. Matlab publish yazıların üretiminde yazı rengi özelliğini kontrol etmeme izin vermiyor. Yazılardaki problem kodlar sadece çok açık rekli. Eğer yazıları fare ile seçerek okuyrsanız yazılar okunabilir oluyor.