showwav.m

some simulation on digital FM,as well as some dsp implemention

function [ output_args ] = showwav( input_args )
%UNTITLED1 Summary of this function goes here
%   Detailed explanation goes here

% 原始信号
[data,fs,nbits]=wavread('recycle.wav');
N=length(data);
fs=38000;
datafft=fftshift(abs(fft(data)));
fx=-0.5:1/N:0.5-1/N;
fx=fx.*fs/1000;
subplot(5,2,1);
plot(fx,datafft);
title('原始信号，38k 采样率');
subplot(5,2,2);
plot(data);
title('时域波形');

% 0 值插值
datainsert(4*N)=0;
k=1;
while k<=N
    datainsert(4*(k-1)+1)=data(k);
    k=k+1;
end
datafft=fftshift(abs(fft(datainsert)));
fx2=-0.5:1/(4*N):0.5-1/(4*N);
fx2=fx2.*(4*fs/1000);
subplot(5,2,3);
plot(fx2,datafft);
title('4 倍插值，采样率变为 38k*4');
subplot(5,2,4);
plot(datainsert);
title('时域波形');

% 低通滤波
b=fir1(64,0.2); % 本来应该是 0.25 的
datafilter=filter(b,1,datainsert);
datafft=fftshift(abs(fft(datafilter)));
subplot(5,2,5);
plot(fx2,datafft);
title('低通滤波后');
subplot(5,2,6);
plot(datafilter);
title('时域波形');

% 调制
w=2*pi*fs/(4*fs);
s=sin(w*[0:4*N-1]);
datas=datafilter.*s;
datafft=fftshift(abs(fft(datas)));
subplot(5,2,7);
plot(fx2,datafft);
title('调幅调制到 38khz');
subplot(5,2,8);
plot(datas);
title('时域波形');

% 插入导频
w=2*pi*(fs/2)/(4*fs);
plevel=0.0003;
s=plevel*sin(w*[0:4*N-1]);
datasum=datafilter+datas+s;
datafft=fftshift(abs(fft(datasum)));
subplot(5,2,9);
plot(fx2,datafft);
title('插入导频和左+右信号');
subplot(5,2,10);
plot(datasum);
title('时域波形');