n=100000
fs=1000
t=(0:n-1)'/fs;
y=cos(2*pi*60*t);
ya=(1+0.1*cos(2*pi*22.2*t));
ch1=y.*ya;
figure
plot (t,y,t,ya)
plot (t,y,t,ch1)
nfft=10000;
nf=nfft/2 + 1;
noverlap=nfft/2;
[P11 F]=pwelch(ch1,hanning(nfft),noverlap,nfft,fs,'onesided');
figure
loglog (F(2:nf),sqrt(P11(2:nf)))
xlabel 'Frequency  (Hz)'
ylabel 'Amplitude  (V/sqrt(Hz))'
