用MATLAB对MIMO MFSK系统进行仿真，画出不同参数下信噪比随误码率的变化曲线，其中参入如下：采样率为20MHz,带宽为2.5MHz，符号速率为625KHz,调制阶数为4，收发天线数目分别为四根发射天线和四根接收天线，莱斯因子分别取5和10，多普勒频移分别取10Hz,4Hz,200Hz,一帧符号数目为1020。

以下是一个简单的MATLAB仿真代码，用于对MIMO MFSK系统进行仿真并绘制误码率随信噪比变化的曲线。代码中的参数可以根据题目中的要求进行调整：

% 信道参数
numTx = 4; % 发送天线数目
numRx = 4; % 接收天线数目
rayleighFactor = [5, 10]; % 莱斯因子
dopplerFreq = [10, 4, 200]; % 多普勒频移

% 信号参数
fs = 20e6; % 采样率
fc = 0; % 中心频率
bw = 2.5e6; % 带宽
symbolRate = 625e3; % 符号速率
numSymbols = 1020; % 一帧符号数目
modOrder = 4; % 调制阶数
modFreq = 1/numSymbols; % 调制频率

% 生成MFSK信号
t = linspace(0, numSymbols/symbolRate, numSymbols*fs/symbolRate);
f = linspace(-bw/2, bw/2, modOrder);
symbols = randi([0, modOrder-1], [numTx, numSymbols]);
modulated = zeros(numTx, length(t));
for i = 1:numTx
    for j = 1:numSymbols
        idx = symbols(i, j) + 1;
        modulated(i, (j-1)*fs/symbolRate+1:j*fs/symbolRate) = cos(2*pi*(f(idx)+fc)*t((j-1)*fs/symbolRate+1:j*fs/symbolRate));
    end
end

% 添加噪声和多径衰落
EbNo = linspace(0, 30, 31); % 信噪比范围
ber = zeros(length(rayleighFactor), length(dopplerFreq), length(EbNo));
for k = 1:length(rayleighFactor)
    for l = 1:length(dopplerFreq)
        for m = 1:length(EbNo)
            % 计算噪声功率
            snr = 10^(EbNo(m)/10);
            noisePower = bw*fs*1/symbolRate*10^(-snr/10);
            
            % 生成多径信道
            channel = rayleighchan(1/symbolRate, dopplerFreq(l), [0 1/fs], [0 -10*log10(rayleighFactor(k))]);
            channel.StoreHistory = true;
            fading = filter(channel, randn(numTx, numSymbols));
            
            % 添加噪声和多径衰落
            received = fading.*modulated + sqrt(noisePower/2)*(randn(numRx, length(t)) + 1i*randn(numRx, length(t)));
            
            % 进行信道估计
            H = zeros(numTx, numRx);
            for i = 1:numSymbols
                H = H + received(:, (i-1)*fs/symbolRate+1:i*fs/symbolRate)*modulated(:, (i-1)*fs/symbolRate+1:i*fs/symbolRate)'/(numSymbols*fs/symbolRate);
            end
            H = H/numSymbols;
            
            % 进行信号检测
            detected = zeros(numRx, numSymbols);
            for i = 1:numSymbols
                detected(:, i) = pinv(H)*received(:, (i-1)*fs/symbolRate+1:i*fs/symbolRate);
            end
            
            % 计算误码率
            [~, errors] = biterr(symbols(:), modQAM(:));
            ber(k, l, m) = errors/length(symbols(:));
        end
    end
end

% 绘制误码率随信噪比变化的曲线
figure;
for k = 1:length(rayleighFactor)
    for l = 1:length(dopplerFreq)
        semilogy(EbNo, squeeze(ber(k, l, :)), 'o-', 'DisplayName', sprintf('Rice Factor=%d, Doppler Shift=%d', rayleighFactor(k), dopplerFreq(l)));
        hold on;
    end
end
grid on;
xlabel('Eb/No (dB)');
ylabel('BER');
legend('show');

这段代码中使用了MATLAB内置的rayleighchan函数生成莱斯信道，使用了pinv函数进行逆矩阵计算进行信号检测，使用了biterr函数计算误码率。由于本问题涉及到较多的参数和复杂的计算，具体的代码实现可能需要进行优化和简化。