pfb

 Example 5: Practical spectrum analysis of a sinusoidal signal

0001 % Example 5: Practical spectrum analysis of a sinusoidal signal
0002 
0003 % Analysis parameters:
0004 M = 31;         % Window length (we'll use a "Hanning window")
0005 N = 64;         % FFT length (zero padding around a factor of 2)
0006 M = 8192;         % Window length (we'll use a "Hanning window")
0007 N = 8192;         % FFT length (zero padding around a factor of 2)
0008 
0009 
0010 % Signal parameters:
0011 wxT = 2*pi/4;   % Sinusoid frequency in rad/sample (1/4 sampling rate)
0012 A = 1;          % Sinusoid amplitude
0013 phix = 0;       % Sinusoid phase
0014 
0015 % Compute the signal x:
0016 n = [0:N-1];    % time indices for sinusoid and FFT
0017 x = A * exp(j*wxT*n+phix); % the complex sinusoid itself: [1,j,-1,-j,1,j,...]
0018 
0019 % Compute Hanning window:
0020 nm = [0:M-1];   % time indices for window computation
0021 %w = (1/M) * (cos((pi/M)*(nm-(M-1)/2))).^2;  % Hanning window = "raised cosine"
0022 w=hamming(M)'
0023 w=blackman(M)'
0024 w=w.*sinc([-M/2:M/2])
0025 
0026 
0027 fs = 500e6;
0028 cutoff =8e6;
0029 t = -256:256;  % This will be a 513-tap filter
0030 r = 2*cutoff/fs;
0031 pfbTaps=2;
0032 unPFBFFTsize=128;
0033 PFBFFTSize=128.*pfbTaps;
0034 t=[-8192./2:8192./2-1];
0035 t=[-unPFBFFTsize./2:unPFBFFTsize./2-1];
0036 tpfb=[-PFBFFTSize./2:PFBFFTSize./2-1];
0037 
0038 wxT = 200e6;   % Sinusoid frequency in rad/sample (1/4 sampling rate)
0039 A = 1;          % Sinusoid amplitude
0040 phix = 0;       % Sinusoid phase
0041 x = A * exp(j*wxT*t+phix); % the complex sinusoid itself: [1,j,-1,-j,1,j,...]
0042 xpfb = A * exp(j*wxT*tpfb+phix); % the complex sinusoid itself: [1,j,-1,-j,1,j,...]
0043 
0044 blackmanfft=fftshift(fft([blackman(length(t))' zeros(1,65536)]));
0045 hammingfft=fftshift(fft([hamming(length(t))' zeros(1,65536)]));
0046 Bblackman=[ x.*sinc(r*t).*r .* blackman(length(t))',zeros(1,65536)];
0047 Bsquare = [x.*sinc(r*t).*r ,zeros(1,65536)];
0048 Bfft=[x,zeros(1,65536)];
0049 Bhamming= [x.*sinc(r*t).*r .* hamming(length(t))' zeros(1,65536)];
0050 Bhammingpfb= [xpfb.*sinc(r*tpfb).*r .* hamming(length(tpfb))'];
0051 
0052 
0053 v=zeros(pfbTaps,length(t))
0054 for i=1:pfbTaps
0055     v(i,:)=Bhammingpfb(1+(i-1).*length(t):(i).*length(t))
0056 end
0057 
0058 BhammingpfbAdded=[v(1,:)+v(2,:) zeros(1,65536)];
0059 
0060 
0061 fftBblackman=fftshift(fft(Bblackman));
0062 fftBsquare=fftshift(fft(Bsquare));
0063 fftBshamming=fftshift(fft(Bhamming));
0064 fftBpfbhamming=fftshift(fft(BhammingpfbAdded));
0065 fftB=fftshift(fft(Bfft,128));
0066 fsS=linspace(1,fs,length(fftBblackman));
0067 import measuredpfbResponse.mat
0068 
0069 
0070 
0071 
0072 fsSquare=linspace(1,fs,length(fftB))
0073 
0074 plot(fsSquare,20*log10(abs(fftB)./max(abs(fftB))))
0075 
0076 plot(fsS,20*log10(abs(fftBshamming)))
0077 hold all
0078 plot(fsS,20*log10(abs(fftBsquare)))
0079 plot(fsS,20*log10(abs(fftBpfbhamming)))
0080 
0081 
0082 wzp = [w,zeros(1,N-M)]; % zero-pad out to the length of x
0083 xw = x .* wzp;          % apply the window w to the signal x
0084 
0085 % Display real part of windowed signal and the Hanning window:
0086 plot(n,wzp,'-'); hold on;
0087 plot(n,real(xw),'*'); 
0088 title('Hanning Window and Windowed, Zero-Padded, Sinusoid (Real Part)'); 
0089 xlabel('Time (samples)'); ylabel('Amplitude'); hold off;
0090 %disp 'pausing for RETURN (check the plot). . .'; pause
0091 %print -deps eps/hanning.eps;
0092 
0093 
0094 
0095 %plot(fn,abs(X),'*'); hold on; plot(fni,abs(Xwi));
0096 % Compute the spectrum and its various alternative forms
0097 %
0098 Xw = fft(xw);              % FFT of windowed data
0099 fn = [0:1.0/N:1-1.0/N];    % Normalized frequency axis
0100 spec = 20*log10(abs(Xw));  % Spectral magnitude in dB
0101 % Since the zeros go to minus infinity, clip at -100 dB:
0102 spec = max(spec,-100*ones(1,length(spec)));
0103 phs = angle(Xw);           % Spectral phase in radians
0104 phsu = unwrap(phs);        % Unwrapped spectral phase (using matlab function)
0105 
0106 Nzp = 16;                   % Zero-padding factor
0107 Nfft = N*Nzp;               % Increased FFT size
0108 xwi = [xw,zeros(1,Nfft-N)]; % New zero-padded FFT buffer
0109 Xwi = fft(xwi);             % Take the FFT
0110 fni = [0:1.0/Nfft:1.0-1.0/Nfft]; % Normalized frequency axis
0111 speci = 20*log10(abs(Xwi));      % Interpolated spectral magnitude in dB
0112 speci = max(speci,-100*ones(1,length(speci))); % clip at -100 dB
0113 phsi = angle(Xwi);          % Phase
0114 phsiu = unwrap(phsi);       % Unwrapped phase
0115 
0116 xi = [x,zeros(1,Nfft-N)];
0117 Xi = fft(xi);  
0118 specxi = 20*log10(abs(Xi));      % Interpolated spectral magnitude in dB
0119 specxi = max(specxi,-20*ones(1,length(specxi))); % clip at -100 dB
0120 phsi = angle(Xi);          % Phase
0121 phsiu = unwrap(phsi);       % Unwrapped phase
0122 
0123 
0124 
0125 
0126 % Plot spectral magnitude
0127 %
0128 
0129 plot(fni,abs(Xwi));
0130 title('Spectral Magnitude'); 
0131 xlabel('Normalized Frequency (cycles per sample))'); 
0132 ylabel('Amplitude (Linear)');
0133 
0134 
0135 plot(fni,abs(Xi));
0136 hold all
0137 plot(fni,abs(Xwi));
0138 title('Spectral Magnitude'); 
0139 xlabel('Normalized Frequency (cycles per sample))'); 
0140 ylabel('Amplitude (Linear)');
0141 
0142 plot(fn,spec,'*'); hold all; 
0143 
0144 fni=[1:length(abs(Xi))]./Nzp;
0145 plot(fni,20*log10(abs(Xi)./max(abs(Xi))));
0146 hold all
0147 plot(fni,20*log10(abs(Xwi)./max(abs(Xwi))));
0148 title('Spectral Magnitude (dB)'); 
0149 xlabel('Normalized Frequency (cycles per sample))'); 
0150 ylabel('Magnitude (dB)');
0151 xlim([2040 2056])
0152 ylim([-60 0])
0153 
0154 
0155 %%%Now we do it with pfb
0156 
0157 pfbTaps=4;
0158 
0159 tap1=xw(1:M./pfbTaps);
0160 tap2=xw(M./pfbTaps+1:2.*M./pfbTaps);
0161 tap3=xw(2.*M./pfbTaps+1:3.*M./pfbTaps);
0162 tap4=xw(3.*M./pfbTaps+1:4.*M./pfbTaps);
0163 xwPfb=tap1+tap2+tap3+tap4;
0164 
0165 fftSize=2048
0166 pfbTaps=2
0167 dataWindowSize=fftSize*pfbTaps
0168 M=dataWindowSize
0169 tap1=xw(1:M./pfbTaps);
0170 tap2=xw(M./pfbTaps+1:2.*M./pfbTaps);
0171 xwPfb2=tap1+tap2;
0172 
0173 Npoly=8192 %size of the polyphase data
0174 Nzp = 16;                   % Zero-padding factor
0175 Nfft = Npoly*Nzp;               % Increased FFT size
0176 xwi = [xw,zeros(1,Nfft-Npoly)]; % New zero-padded FFT buffer
0177 Xwi = fft(xwi);             % Take the FFT
0178 fnixw= [0:1.0/Nfft:1.0-1.0/Nfft]; % Normalized frequency axis
0179 
0180 XwiDec=Xwi(1:pfbTaps:length(Xwi))
0181 fnixw= [0:1.0/Nfft:1.0-1.0/Nfft]; % Normalized frequency axis
0182 
0183 
0184 Npoly=8192 %size of the polyphase data
0185 Nzp = 16;                   % Zero-padding factor
0186 Nfft = Npoly*Nzp;               % Increased FFT size
0187 x = [x,zeros(1,Nfft-Npoly)]; % New zero-padded FFT buffer
0188 X = fft(x);             % Take the FFT
0189 fnix= [0:1.0/Nfft:1.0-1.0/Nfft]; % Normalized frequency axis
0190 
0191 Npoly=2048 %size of the polyphase data
0192 XwPfB = fft(xwPfb);  
0193 Nzp = 16;                   % Zero-padding factor
0194 Nfft = Npoly*Nzp;               % Increased FFT size
0195 xwpfbi = [xwPfb,zeros(1,Nfft-Npoly)]; % New zero-padded FFT buffer
0196 XwPfbi = fft(xwpfbi);             % Take the FFT
0197 fniXwp = [0:1.0/Nfft:1.0-1.0/Nfft]; % Normalized frequency axis
0198 
0199 Npoly=2048 %size of the polyphase data
0200 XwPfB2 = fft(xwPfb2);  
0201 Nzp = 16;                   % Zero-padding factor
0202 Nfft = Npoly*Nzp;               % Increased FFT size
0203 xwpfb2i = [xwPfb2,zeros(1,Nfft-Npoly)]; % New zero-padded FFT buffer
0204 XwPfb2i = fft(xwpfb2i);             % Take the FFT
0205 fni = [0:1.0/Nfft:1.0-1.0/Nfft]; % Normalized frequency axis
0206 
0207 
0208 
0209 Nfft=32768
0210 xwi = [xw,zeros(1,(Nfft-length(xw)))]; % New zero-padded FFT buffer
0211 Xwi = fft(xwi);             % Take the FFT
0212 
0213 
0214 
0215 fni=[1:length(abs(XwPfbi))]./Nzp;
0216 
0217 
0218 plot(fni,20*log10(abs(XwPfbi)./max(abs(XwPfbi))));
0219 hold all
0220 plot(fni,20*log10(abs(XwPfb2i)./max(abs(XwPfb2i))));
0221 
0222 plot(fnixw,20*log10(abs(Xwi)./max(abs(Xwi))));
0223 hold all
0224 plot(fnix,20*log10(abs(X)./max(abs(X))));
0225 
0226 
0227 
0228 
0229 
0230 
0231 plot(20*log10(abs(XwPfB)./max(abs(XwPfB))),'r');
0232 xlim([500 520])
0233 
0234 plot(20*log10(abs(XwPfB)))
0235 
0236 
0237 xwPfbi=resample(xwPfb,16,1);
0238 XwPfbi = fft(xwpfbi);             % Take the FFT
0239 
0240 plot(fni(1:4:length(fni)),20*log10(abs(Xwi(1:4:length(fni))./max(abs(Xwi)))));
0241 
0242