[y,Fs,bits]=wavread('DSP(1213)_A.wav'); % read sound file
N=round(5.56/0.001); % set the order of the filter
%F1=0.07;
F2=0.29; %ÓÉÐźÅDSP(1213)_C.wavµÄƵÂʳɷֶ¨500-2000%wp=fp/(Fs/2);
Fs=14000
%F1=0.10;F2=0.38; %ÊÊÓÚDSP(1213)_B.wav ÒòΪ¸ÃÐźŵÄƵÂÊ´óÖÂΪ600¡ª¡ª2300 Fs=12000
F1=0.12;F2=0.33;%ÊÊÓÚDSP(1213)_A.wav
ÒòΪ¸ÃÐźŵÄƵÂÊ´óÖÂΪ700¡ª¡ª1900/19 Fs=11250
b=fir1(N,[F1,F2],blackman(N+1)); % using blackman window to implement FIR filter
%b=fir1(N,[F1,F2],hanning(N+1));
Y=filter(b,1,y); % output te signal after filtering
%freqz(b,1,512,Fs);
figure
freqz(b,1);%normalized frequency
title('the filter response(magnitude and Phase)without using scaling factor after using blackman window');
%title('the filter response(magnitude and Phase)without using scaling factor');%after using blackman window');
t=(0:length(y)-1)/Fs;
figure;
subplot(2,1,1);
plot(t,y); % plot original signal in time domain
title('Input signal in time domain');
xlabel('Time(second)');
ylabel('Magnitude');
subplot(2,1,2);
plot(t,Y); % plot output signal in time domain
title('Output signal in time domain');
xlabel('Time(second)');
ylabel('Magnitude');
figure;
fft_y=fft(y);
fft_Y=fft(Y);
f=Fs*(0:length(y)-1)/length(y);
subplot(2,1,1);
plot(f,20*log10(abs(fft_y(1:length(y))))); % plot original signal in frequency domain
%plot(f,abs(fft_y(1:length(y))));
title('Input signal in frequency domain');
xlabel('Frequency(Hz)');
ylabel('Amplitude(dB)');
subplot(2,1,2);
plot(f,20*log10(abs(fft_Y(1:length(y))))); % plot output signal in frequency domain
%plot(f,abs(fft_Y(1:length(y))));
title('Output signal in frequency domain');
xlabel('Frequency(Hz)');
ylabel('Amplitude(dB)');
figure;
subplot(2,1,1);
%specgram(y,[],Fs); %plot original signal in mixed domain
spectrogram(y,[],Fs);
colorbar;
title('Input signal energy spectrum');
%xlabel('Time(second)');
%ylabel('Frequency(Hz)');
grid on;
axis tight;
subplot(2,1,2);
%specgram(Y,[],Fs); %plot output signal in mixed domain
spectrogram(Y,[],Fs);
colorbar;
title('Output signal energy spectrum');
%xlabel('Time(second)');
%ylabel('Frequency(Hz)');
grid on;
axis tight;
factor=(2^15-1)/max(b);
m=round(log2(factor)); %decide m
scalf=2^m; %decide scaling factor based on m
new_b1=round(b*scalf);
new_b2=new_b1/scalf; %new coefficient
figure;
freqz(new_b2,1); %plot new filter response
%freqz(new_b2,1,512,Fs);
%title('the filter response(magnitude and Phase)with using scaling factor') figure;
Y2=filter(new_b2,1,y);
fft_Y=fft(Y,2048);
fft_Y2=fft(Y2,2048);
f=Fs*(0:1023)/2048;
plot(f,20*log10(abs(fft_Y(1:1024))),'r');% plot frequency responce without scalling factor
hold on;
plot(f,20*log10(abs(fft_Y2(1:1024)))); %plot frequency responce with scalling factor
title('frequency responce with and without scaling factor');
xlabel('Frequency(Hz)');
ylabel('Amplitude(dB)');
legend('without scaling factor','with scaling factor'); %scalling factor Ëõ·Å±ÈÀýÒò×Ó
wavplay(Y.*200,Fs); %make sound
wavwrite(Y.*200,Fs,'final sound.wav'); %stroe the nex sound file
wavplay(Y2.*200,Fs);
wavwrite(Y2.*200,Fs,'final sound2.wav');。