[問題] 訊號處理程式問題 求救已回收
clc;
clear;
close all
load('D:\課業相關\XXXXXX\程式\matlab.mat')
fs=500;
%資料Sheet1檔名
N=length(Sheet1)
t=(0:1/fs:(length(Sheet1)/fs-1/fs))'
T=1/fs;
m=0:N-1;
fd=1/(N*T);
fy=m.*fd;
%去除直流和多項是飄移
pl_Signal=detrend(Sheet1)
%原始圖
figure(1)
plot(t,Sheet1)
title('原始圖')
xlabel('時間');
ylabel('電壓值');
%原始圖FFT
a=fft(Sheet1)
figure(2)
plot(fy,abs(a))
title('頻譜圖')
xlabel('頻率');
ylabel('強度值');
%去除直流和多項是飄移圖
figure(3)
plot(t,pl_Signal)
title('多項式飄移圖')
xlabel('時間');
ylabel('電壓');
%去除直流和多項是飄移圖FFT
a=fft(pl_Signal)
figure(4)
plot(fy,abs(a))
title('多項式飄移頻譜圖')
xlabel('頻率');
ylabel('強度值');
%加窗
add_win=hamming(N);
y_win=pl_Signal.*(add_win)
figure(5)
plot(t,y_win)
title('hamming加窗圖')
xlabel('時間');
ylabel('電壓值');
%加窗FFT
a=fft(y_win)
figure(6)
plot(fy,abs(a))
title('hamming加窗頻譜圖')
xlabel('頻率');
ylabel('強度值');
%高通 1階1HZ
[b,a]=butter(1,0.1/(fs/2),'high');
high=filter(b,a,y_win);
figure(7)
plot(t,high)
title('高通')
xlabel('時間');
ylabel('電壓值');
figure(8)
a=fft(high)
plot(fy,abs(a))
title('高通後頻譜圖')
xlabel('頻率');
ylabel('強度值');
%低通 5階40HZ
[b,a]=butter(5,40/(fs/2),'low');
low=filter(b,a, high);
figure(9)
plot(t,low)
title('低通')
xlabel('時間');
ylabel('電壓值');
figure(10)
a=fft(low)
plot(fy,abs(a))
title('低通')
xlabel('頻率');
ylabel('強度值');
%10點平均FIR濾波器
fs=500;
t=(0:1/fs:(length(Sheet1)/fs-1/fs))';
N=length(t);
T=1/fs;
fd=1/(N*T);
m=0:N-1;
fy=m.*fd;
x=pl_Signal
b=ones(1,10)/10; %10點平均的FIR濾波器
y1=filtfilt(b,1,x);
y2=filter(b,1,x);
figure(11)
plot(t,y1,'r')
xlabel('Time (sec)')
ylabel('Amplitude')
title('Anti-Causal, zero-phase filter implementation')
figure(12)
plot(t,y2,'g')
xlabel('Time (sec)')
ylabel('Amplitude')
title('Anti-Causal, zero-phase filter implementation')
%Chebyshev 型I
[n,Wn]=cheb1ord(40/(fs/2),50/(fs/2),1,30);
[b,a]=cheby1(n,1,Wn);
figure(13)
freqz(b,a,512,fs)
r=filter(b,a,y_win)
figure(14)
plot(t,r)
r1=fft(r)
figure(15)
plot(fy,abs(r1))
%Chebyshev 型II
[n,Wn]=cheb2ord(40/(fs/2),50/(fs/2),1,30);
[b,a]=cheby2(n,1,Wn);
figure(16)
freqz(b,,512,fs)
r=filter(b,a,y_win)
figure(17)
plot(t,r)
r1=fft(r)
figure(18)
plot(fy,abs(r1))
%Elliptic
[n,Wn]=ellipord(40/(fs/2),50/(fs/2),1,30);
[b,a]=ellip(n,1,60,Wn);
figure(19)
freqz(b,a,512,fs)
r=filter(b,a,y_win)
figure(20)
plot(t,r)
r1=fft(r)
figure(21)
plot(fy,abs(r1))
MATLAB說
Error in ==> W1 at 121
[n,Wn]=cheb1ord(40/(fs/2),50/(fs/2),1,30);
問題出在哪了??
快瘋掉~~希望各位幫幫忙了
感謝!!
--
人生總會遇到畜生?!
--
※ 發信站: 批踢踢實業坊(ptt.cc)
◆ From: 122.117.189.65
→
04/30 11:11, , 1F
04/30 11:11, 1F
MATLAB 近期熱門文章
PTT數位生活區 即時熱門文章
