close all;

fs = 1e6;           % sampling frequency
cycles = 43.4;      % number of full cycles in sinusoid
N = 2^10;           % number of time steps in analysis
                    % (power of 2 speeeds up analyis)
fx = fs*cycles/N    % signal frequency

t = linspace(0, (N-1)/fs, N);
y = sin(2*pi*fx*t);
w = nuttallwin(N)';
w = w/sqrt(var(w)*2);
s = abs(fft(y)/N*2);
s = s(1:N/2);
f = (0:length(s)-1) / N;
sw= abs(fft(y .* w)/N*2);
sw= sw(1:N/2);

th = 1e-4;
amp_with_window = sw(sw>th)
number_of_bins = sum(sw>th)
power_in_these_bins = sum(sw(sw>th) .^ 2)

figure(1);
plot(w);
axis([ 1 N min(w) 1.1*max(w) ]);
fixfig;

figure(2);
ww= abs(fft(w)/mean(w)/N);
ww= 20 * log10(ww(1:20));
plot(1:length(ww), ww, 'rx');
xlabel('DFT Bin');
ylabel('Normalized Aplitude');
axis([1 length(ww) min(ww) max(ww) ]);
fixfig;

figure(3);
subplot(2, 1, 1);
plot(f, 20*log10(s));
xlabel('Frequency  [ f / f_s ]');
ylabel('Spectrum not Windowed  [ dBFS ]');
fixfig;
subplot(2, 1, 2);
plot(f, 20*log10(sw));
xlabel('Frequency  [ f / f_s ]');
ylabel('Windowed Spectrum  [ dBFS ]');
fixfig;

figure(4);
subplot(2, 1, 1);
plot(t, y);
xlabel('Time');
ylabel('Signal Amplitude');
axis([ min(t) max(t) 1.1*min(y) 1.1*max(y) ]);
subplot(2, 1, 2);
plot(t, y .* w);
xlabel('Time');
ylabel('Windowed Signal Amplitude');
axis([ min(t) max(t) 1.1*min(y .* w) 1.1*max(y .* w) ]);