# FFT of quantized signal

35 vues (au cours des 30 derniers jours)
Maruthi prasanna C le 2 Fév 2020
Commenté : Addy le 17 Août 2022
I am having difficulty in plotting FFT of quantized signal of an ideal qunatizer in MATLAB . When I plot the FFT in dB Scale , I am getting infinity as the value & figure is blank . Please help me, I have posted the Code below
fs = 1e6;
fx = 50e3;
Afs = 1;
N = 2^11;
% time vector
t = linspace(0, (N-1)/fs, N);
% input signal
y = Afs * cos(2*pi*fx*t);
% spectrum
s = 20 * log10(abs(fft(y)/N/Afs*2));
% drop redundant half
s = s(1:N/2);
% frequency vector (normalized to fs)
f = (0:length(s)-1);
figure(1)
plot(y);
title('input signal ')
ylabel('y');
xlabel('frequency bins');
B = 2; % no of bits
delta = 2/2^B;
fs = 1e6;
% Number of cycles in test
cycles = 67;
%Make N/cycles non-integer! accomplished by choosing cyclesÆprime #
%N=power of 2 speeds up analysis
% N = 2^10;
%signal frequency
fx = fs*cycles/N
y = Afs * cos(2*pi*fx*t);
s = 20 * log10(abs(fft(y)/N/Afs*2));
f = (0:length(s)-1);
figure(2)
plot(f,s);
title('spectrum of input signal ');
ylabel('y')
xlabel('frequency bins')
y1=round(y/delta)*delta;
figure(3)
plot(y,'b');
hold on
stem(y1,'r')
hold on
plot(y-y1,'g')
title(sprintf('Signal, Quantized signal and Error for %g bits, %g quantization levels',B,2^B));
xlabel('frequency');
hold off
y2=20*log10((abs(fft(y1)/N*2)));
figure(10)
f1 = (0:length(y2)-1);
plot(f1 ,y2)
##### 4 commentairesAfficher 3 commentaires plus anciensMasquer 3 commentaires plus anciens
Chunru le 11 Août 2021
In theory, quantization can have any integer number of levels (not necessarily 2^B, though 2^B is prefered). The quantization can be non-uniform as well.

Connectez-vous pour commenter.

### Réponses (2)

Chunru le 11 Août 2021
Quantizing the ideal sinusoidal signal will produce a perodical signal, which may have some weird spectrum with 0 at some frequencies causing the display problem. To study the quatization effect, you can use some other signals, for example chirp, pink noise. The following code add in an offset so that the spectrum is no longer "weird".
fs = 1e6;
fx = 50e3;
Afs = 1;
N = 2^11;
% time vector
t = linspace(0, (N-1)/fs, N);
% input signal
y = Afs * cos(2*pi*fx*t);
% spectrum
s = 20 * log10(abs(fft(y)/N/Afs*2));
% drop redundant half
s = s(1:N/2);
% frequency vector (normalized to fs)
f = (0:length(s)-1);
figure(1)
plot(y);
title('input signal ')
ylabel('y');
xlabel('frequency bins'); B = 2; % no of bits
delta = 2/2^B;
fs = 1e6;
% Number of cycles in test
cycles = 167;
%Make N/cycles non-integer! accomplished by choosing cyclesÆprime #
%N=power of 2 speeds up analysis
% N = 2^10;
%signal frequency
fx = fs*cycles/N
fx = 8.1543e+04
y = Afs * cos(2*pi*fx*t) + 0.123; % add a small offset
s = 20 * log10(abs(fft(y)/N/Afs*2));
f = (0:length(s)-1);
figure(2)
plot(f,s);
title('spectrum of input signal ');
ylabel('y')
xlabel('frequency bins') y1=round(y/delta)*delta; y1'
ans = 2048×1
1.0000 1.0000 0.5000 0 -0.5000 -0.5000 -1.0000 -1.0000 -0.5000 0
figure(3)
plot(y,'b');
hold on
stem(y1,'r')
hold on
plot(y-y1,'g')
title(sprintf('Signal, Quantized signal and Error for %g bits, %g quantization levels',B,2^B));
xlabel('frequency');
hold off y2=20*log10((abs(fft(y1)/N*2)))'
y2 = 2048×1
-13.8447 -58.6345 -62.4346 -74.3996 -103.6032 -59.0645 -60.1435 -59.3283 -45.3150 -38.9921
%y2=((abs(fft(y1)/N*2)))'
figure(10)
f1 = (0:length(y2)-1);
plot(f1 ,y2) ##### 0 commentairesAfficher -1 commentaires plus anciensMasquer -1 commentaires plus anciens

Connectez-vous pour commenter.

Yazan le 11 Août 2021
There are some issues with your code. See my modifications and comments below. You cannot calculate SNR, because you did not add noise to any signal (unless I'm missing something here).
clc, clear, close all
fs = 1e6;
fx = 50e3;
Afs = 1;
N = 2^11;
% time vector
t = linspace(0, (N-1)/fs, N);
% input signal
y = Afs * cos(2*pi*fx*t);
figure('Units', 'normalized', 'Position', [0.05 0.15 0.85 0.7])
subplot(3,2,1)
plot(t, y), title('Signal y'); axis tight
xlabel('Time - sec');
ylabel('Amplitude');
% spectrum
Nfft = N;
s = 1/Nfft * abs(fft(y, Nfft));
f = linspace(0, fs/2, Nfft/2+1);
subplot(3,2,2)
s = s(1:Nfft/2+1).^2;
sy = [s(1), 2*s(2:end-1), s(end)];
plot(f, pow2db(sy)), grid minor
title('One-sided power spectrum of y'), xlabel('Frequency - Hz');
ylabel('dB')
% or you can use the periodogram function with rect window
% you'll get exactly the same result
% [s2, f2] = periodogram(y, rectwin(length(y)), Nfft, fs, 'onesided', 'power');
% hold on, plot(f2, pow2db(s2))
% signal 2
B = 2;
delta = 2/2^B;
fs = 1e6;
cycles = 67;
fx = fs*cycles/N;
yy = Afs * cos(2*pi*fx*t);
subplot(3,2,3)
plot(t, yy), title('Signal yy'); axis tight
xlabel('Time - sec');
ylabel('Amplitude');
% spectrum
Nfft = N;
s = 1/Nfft * abs(fft(yy, Nfft));
s = s(1:Nfft/2+1).^2;
syy = [s(1), 2*s(2:end-1), s(end)];
subplot(3,2,4)
plot(f, pow2db(syy)), grid minor
title('One-sided power spectrum of yy'), xlabel('Frequency - Hz');
ylabel('dB')
% quantized signal
y1 = round(y/delta)*delta;
subplot(3,2,5)
n = 100;
plot(t(1:n), y1(1:n)), title('Signal y1'); axis tight
xlabel('Time - sec');
ylabel('Amplitude'); grid minor
hold on,
plot(t(1:n), y(1:n), '--')
err = y-y1;
plot(t(1:n), err(1:n))
legend({'Quantized', 'Original', 'Error'})
s = 1/Nfft * abs(fft(y1, Nfft));
s = s(1:Nfft/2+1).^2;
sy1 = [s(1), 2*s(2:end-1), s(end)];
subplot(3,2,6)
plot(f, pow2db(sy1)), grid minor
title('One-sided power spectrum of y1'), xlabel('Frequency - Hz');
ylabel('dB') % estimate power
pRMSy = mag2db(rms(y));
pRMSyy = mag2db(rms(yy));
pRMSy1 = mag2db(rms(y1));
fprintf('Power of signal y = %g dBW\n', pRMSy)
Power of signal y = -3.01267 dBW
fprintf('Power of signal yy = %g dBW\n', pRMSyy)
Power of signal yy = -3.0103 dBW
fprintf('Power of signal y1 = %g dBW\n', pRMSy1)
Power of signal y1 = -2.22132 dBW
##### 1 commentaireAfficher AucuneMasquer Aucune
Hey in your code, 2 bit quantization, there are 5 levels.
you can see when you do "unique(y1)"
But for 2 bit quantization, there is only 4 levels. 00, 01, 10, 11.
So, dividing your signal into 4 parts is necessary (At least that is what I am understanding. Please correct me if I am wrong. I am still learning.
I was working on that 4 level quantization.
clear;clc
A = 4;
f = 10e3;
% fs = 25e6; ts = 1/fs;
ts = 0.2/10/f(end);
t_end = 1/f;
t = 0:ts:t_end;
x = A*sin(2*pi*f*t)';
% x = awgn(x,100);
t = t*1e6;
%%
q_bits = 2;
q_num_levels = 2^q_bits;
q_levels = linspace(min(x),max(x),q_num_levels);
bins = interp1(1:numel(q_levels),q_levels,0.5 + (1:numel(q_levels)-1));
y = discretize(x,[-Inf bins Inf]);
y = normalize(y,'range',[min(q_levels),max(q_levels)]);
y = y-mean(y);
u = unique(y);
%%
bit_levels = 0:q_num_levels-1;
bit_levels = dec2bin(bit_levels,q_bits);
q = diff(q_levels); q = q(1);
q_error = (max(x)-min(x))/(2*q_num_levels);
q_noise_pwr = q/sqrt(12);
%%
figure(1);clf
hold on;box on
set(gcf,'color','k')
set(gca,'color','k')
set(gca, 'YGrid', 'on', 'XGrid', 'off')
set(gca,'GridColor','w')
set(gca,'XColor','w')
set(gca,'YColor','w')
set(gca,'GridAlpha',0.5)
plot(t,x,'c')
plot(t,y,'m');
plot(t,x-y,'-.g')
yticks(q_levels);
yticklabels(bit_levels);
yticklabels([]);
ylim([min(x)+min(x)/10,max(x)+max(x)/10])
xlim([min(t),max(t)])
tit = title(['Example with ',num2str(q_bits),' Bits[',num2str(q_num_levels),' levels] Linear Quantization']);
lgnd = legend( {'\color{cyan} Analog signal',...
'\color{magenta} Quantized signal',...
['\color{green} Quantized Error signal.',newline,...
' Max.Quant.Error = ',num2str(q_error)]});
set(lgnd, 'Color', 'k');
set(tit, 'Color', 'w');
ylabel('Quantization levels')
xlabel('Time - [µs]') Connectez-vous pour commenter.

### Catégories

En savoir plus sur Spectral Measurements dans Help Center et File Exchange

### Community Treasure Hunt

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

Start Hunting!