diff --git a/Project.toml b/Project.toml
index aa10c3f..a8d4ada 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "JudiLing"
 uuid = "b43a184b-0e9d-488b-813a-80fd5dbc9fd8"
 authors = ["Xuefeng Luo", "Maria Heitmeier"]
-version = "0.8.4"
+version = "0.9.0"
 
 [deps]
 BSON = "fbb218c0-5317-5bc6-957e-2ee96dd4b1f0"
diff --git a/README.md b/README.md
index 858a6a3..a476a28 100644
--- a/README.md
+++ b/README.md
@@ -1,8 +1,8 @@
 # JudiLing
 
-[![Stable](https://img.shields.io/badge/docs-stable-blue.svg)](https://MegamindHenry.github.io/JudiLing.jl/stable)
-[![Dev](https://img.shields.io/badge/docs-dev-blue.svg)](https://MegamindHenry.github.io/JudiLing.jl/dev)
-[![Build Status](https://github.com/MegamindHenry/JudiLing.jl/workflows/CI/badge.svg)](https://github.com/MegamindHenry/JudiLing.jl/actions)
+[![Stable](https://img.shields.io/badge/docs-stable-blue.svg)](https://quantling.github.io/JudiLing.jl/stable)
+[![Dev](https://img.shields.io/badge/docs-dev-blue.svg)](https://quantling.github.io/JudiLing.jl/dev)
+[![Build Status](https://github.com/quantling/JudiLing.jl/workflows/CI/badge.svg)](https://github.com/quantling/JudiLing.jl/actions)
 [![codecov](https://codecov.io/gh/MegamindHenry/JudiLing.jl/branch/master/graph/badge.svg)](https://codecov.io/gh/MegamindHenry/JudiLing.jl)
 
 JudiLing: An implementation for Linear Discriminative Learning in Julia
diff --git a/docs/src/index.md b/docs/src/index.md
index 8852d4b..9737e94 100644
--- a/docs/src/index.md
+++ b/docs/src/index.md
@@ -14,11 +14,11 @@ Pkg.add("JudiLing")
 ```
 For brave adventurers, install test version of JudiLing by:
 ```
-julia> Pkg.add(url="https://github.com/MegamindHenry/JudiLing.jl.git")
+julia> Pkg.add(url="https://github.com/quantling/JudiLing.jl.git")
 ```
 Or from the Julia REPL, type `]` to enter the Pkg REPL mode and run
 ```
-pkg> add https://github.com/MegamindHenry/JudiLing.jl.git
+pkg> add https://github.com/quantling/JudiLing.jl.git
 ```
 
 ## Running Julia with multiple threads
diff --git a/docs/src/man/make_cue_matrix.md b/docs/src/man/make_cue_matrix.md
index b6d24ba..e43df5f 100644
--- a/docs/src/man/make_cue_matrix.md
+++ b/docs/src/man/make_cue_matrix.md
@@ -13,5 +13,12 @@ CurrentModule = JudiLing
     make_cue_matrix(data::DataFrame, cue_obj::Cue_Matrix_Struct)
     make_cue_matrix(data_train::DataFrame, data_val::DataFrame)
     make_combined_cue_matrix(data_train, data_val)
+    make_cue_matrix_from_CFBS(features::Vector{Vector{T}};
+                                        pad_val::T = 0.,
+                                        ncol::Union{Missing,Int}=missing) where {T}
+    make_combined_cue_matrix_from_CFBS(features_train::Vector{Vector{T}},
+                                                features_test::Vector{Vector{T}};
+                                                pad_val::T = 0.,
+                                                ncol::Union{Missing,Int}=missing) where {T}
     make_ngrams(tokens, grams, keep_sep, sep_token, start_end_token)
 ```
diff --git a/src/cholesky.jl b/src/cholesky.jl
index 47b2a6f..ec8e056 100644
--- a/src/cholesky.jl
+++ b/src/cholesky.jl
@@ -327,7 +327,7 @@ end
 """
     make_transform_matrix(X::Union{SparseMatrixCSC,Matrix},
                             Y::Union{SparseMatrixCSC,Matrix},
-                            freq::Array{Int64,1})
+                            freq::Union{Array{Int64, 1}, Array{Float64,1}})
 
 Weight X and Y using the frequencies in freq. Then use the Cholesky
 decomposition to calculate the transformation matrix from X to Y,
@@ -336,7 +336,7 @@ where X is a sparse matrix and Y is a sparse matrix.
 # Obligatory Arguments
 - `X::SparseMatrixCSC`: the X matrix, where X is a sparse matrix
 - `Y::SparseMatrixCSC`: the Y matrix, where Y is a sparse matrix
-- `freq::Array{Int64,1}`: list of frequencies of the wordforms in X and Y
+- `freq::Union{Array{Int64, 1}, Array{Float64,1}}`: list of frequencies of the wordforms in X and Y
 
 # Optional Arguments
 - `method::Symbol = :additive`: whether :additive or :multiplicative decomposition is required
diff --git a/src/eval.jl b/src/eval.jl
index 4392266..b01cf7c 100644
--- a/src/eval.jl
+++ b/src/eval.jl
@@ -225,6 +225,7 @@ function accuracy_comprehension(
     Comp_Acc_Struct(dfr, acc, err)
 end
 
+
 """
     eval_SC(SChat::AbstractArray, SC::AbstractArray)
 
@@ -232,6 +233,8 @@ Assess model accuracy on the basis of the correlations of row vectors of Chat an
 C or Shat and S. Ideally the target words have highest correlations on the diagonal
 of the pertinent correlation matrices.
 
+If `freq` is added, token-based accuracy is computed. Token-based accuracy weighs accuracy values according to words' frequency, i.e. if a word has a frequency of 30 and overall there are 3000 tokens (the frequencies of all types sum to 3000), this token's accuracy will contribute 30/3000.
+
 !!! note
     If there are homophones/homographs in the dataset, this evaluation method may be misleading: the predicted vector will be equally correlated with the target vector of both words and the one on the diagonal will not necessarily be selected as the most correlated. In such cases, supplying the dataset and `target_col` is recommended which enables taking into account homophones/homographs.
 
@@ -242,6 +245,7 @@ of the pertinent correlation matrices.
 # Optional Arguments
 - `digits`: the specified number of digits after the decimal place (or before if negative)
 - `R::Bool=false`: if true, pairwise correlation matrix R is return
+- `freq::Union{Missing, Array{Int64, 1}, Array{Float64,1}}=missing`: list of frequencies of the wordforms in X and Y
 
 ```julia
 eval_SC(Chat_train, cue_obj_train.C)
@@ -250,7 +254,8 @@ eval_SC(Shat_train, S_train)
 eval_SC(Shat_val, S_val)
 ```
 """
-function eval_SC(SChat::AbstractArray, SC::AbstractArray; digits=4, R=false)
+function eval_SC(SChat::AbstractArray, SC::AbstractArray; digits=4, R=false,
+    freq::Union{Missing, Array{Int64, 1}, Array{Float64,1}}=missing)
 
     if size(unique(SC, dims=1), 1) != size(SC, 1)
         @warn "eval_SC: The C or S matrix contains duplicate vectors (usually because of homophones/homographs). Supplying the dataset and target column is recommended for a realistic evaluation. See the documentation of this function for more information."
@@ -262,7 +267,12 @@ function eval_SC(SChat::AbstractArray, SC::AbstractArray; digits=4, R=false)
         dims = 2,
     )
     v = [rSC[i[1], i[1]] == rSC[i] ? 1 : 0 for i in argmax(rSC, dims = 2)]
-    acc = round(sum(v) / length(v), digits=digits)
+    if !ismissing(freq)
+        v .*= freq
+        acc = round(sum(v) / sum(freq), digits=digits)
+    else
+        acc = round(sum(v) / length(v), digits=digits)
+    end
     if R
         return acc, rSC
     else
@@ -277,6 +287,8 @@ Assess model accuracy on the basis of the correlations of row vectors of Chat an
 C or Shat and S. Ideally the target words have highest correlations on the diagonal
 of the pertinent correlation matrices.
 
+If `freq` is added, token-based accuracy is computed. Token-based accuracy weighs accuracy values according to words' frequency, i.e. if a word has a frequency of 30 and overall there are 3000 tokens (the frequencies of all types sum to 3000), this token's accuracy will contribute 30/3000.
+
 !!! note
     The order is important. The fist gold standard matrix has to be corresponing
     to the SChat matrix, such as `eval_SC(Shat_train, S_train, S_val)` or `eval_SC(Shat_val, S_val, S_train)`
@@ -292,6 +304,7 @@ of the pertinent correlation matrices.
 # Optional Arguments
 - `digits`: the specified number of digits after the decimal place (or before if negative)
 - `R::Bool=false`: if true, pairwise correlation matrix R is return
+- `freq::Union{Missing, Array{Int64, 1}, Array{Float64,1}}=missing`: list of frequencies of the wordforms in X and Y
 
 ```julia
 eval_SC(Chat_train, cue_obj_train.C, cue_obj_val.C)
@@ -305,12 +318,14 @@ function eval_SC(
     SC::AbstractArray,
     SC_rest::AbstractArray;
     digits = 4,
-    R = false
+    R = false,
+    freq::Union{Missing, Array{Int64, 1}, Array{Float64,1}}=missing
     )
 
-    eval_SC(SChat, vcat(SC, SC_rest); digits=digits, R=R)
+    eval_SC(SChat, vcat(SC, SC_rest); digits=digits, R=R, freq=freq)
 end
 
+
 """
     eval_SC(SChat::AbstractArray, SC::AbstractArray, data::DataFrame, target_col::Union{String, Symbol})
 
@@ -318,6 +333,8 @@ Assess model accuracy on the basis of the correlations of row vectors of Chat an
 C or Shat and S. Ideally the target words have highest correlations on the diagonal
 of the pertinent correlation matrices. Support for homophones.
 
+If `freq` is added, token-based accuracy is computed. Token-based accuracy weighs accuracy values according to words' frequency, i.e. if a word has a frequency of 30 and overall there are 3000 tokens (the frequencies of all types sum to 3000), this token's accuracy will contribute 30/3000.
+
 # Obligatory Arguments
 - `SChat::Union{SparseMatrixCSC, Matrix}`: the Chat or Shat matrix
 - `SC::Union{SparseMatrixCSC, Matrix}`: the C or S matrix
@@ -327,6 +344,7 @@ of the pertinent correlation matrices. Support for homophones.
 # Optional Arguments
 - `digits`: the specified number of digits after the decimal place (or before if negative)
 - `R::Bool=false`: if true, pairwise correlation matrix R is return
+- `freq::Union{Missing, Array{Int64, 1}, Array{Float64,1}}=missing`: list of frequencies of the wordforms in X and Y
 
 ```julia
 eval_SC(Chat_train, cue_obj_train.C, latin, :Word)
@@ -341,7 +359,8 @@ function eval_SC(
     data::DataFrame,
     target_col::Union{String, Symbol};
     digits = 4,
-    R = false
+    R = false,
+    freq::Union{Missing, Array{Int64, 1}, Array{Float64,1}}=missing
     )
 
     rSC = cor(
@@ -353,7 +372,12 @@ function eval_SC(
         data[i[1], target_col] == data[i[2], target_col] ? 1 : 0
         for i in argmax(rSC, dims = 2)
     ]
-    acc = round(sum(v) / length(v), digits=digits)
+    if !ismissing(freq)
+        v .*= freq
+        acc = round(sum(v) / sum(freq), digits=digits)
+    else
+        acc = round(sum(v) / length(v), digits=digits)
+    end
     if R
         return acc, rSC
     else
@@ -368,6 +392,8 @@ Assess model accuracy on the basis of the correlations of row vectors of Chat an
 C or Shat and S. Ideally the target words have highest correlations on the diagonal
 of the pertinent correlation matrices.
 
+If `freq` is added, token-based accuracy is computed. Token-based accuracy weighs accuracy values according to words' frequency, i.e. if a word has a frequency of 30 and overall there are 3000 tokens (the frequencies of all types sum to 3000), this token's accuracy will contribute 30/3000.
+
 !!! note
     The order is important. The fist gold standard matrix has to be corresponing
     to the SChat matrix, such as `eval_SC(Shat_train, S_train, S_val, latin, :Word)`
@@ -384,6 +410,7 @@ of the pertinent correlation matrices.
 # Optional Arguments
 - `digits`: the specified number of digits after the decimal place (or before if negative)
 - `R::Bool=false`: if true, pairwise correlation matrix R is return
+- `freq::Union{Missing, Array{Int64, 1}, Array{Float64,1}}=missing`: list of frequencies of the wordforms in X and Y
 
 ```julia
 eval_SC(Chat_train, cue_obj_train.C, cue_obj_val.C, latin, :Word)
@@ -400,7 +427,8 @@ function eval_SC(
     data_rest::DataFrame,
     target_col::Union{String, Symbol};
     digits = 4,
-    R = false
+    R = false,
+    freq::Union{Missing, Array{Int64, 1}, Array{Float64,1}}=missing
     )
 
     n_data = size(data, 1)
@@ -428,7 +456,8 @@ function eval_SC(
         data_combined,
         target_col,
         digits = digits,
-        R = R
+        R = R,
+        freq=freq
         )
 end
 
diff --git a/src/make_cue_matrix.jl b/src/make_cue_matrix.jl
index e5110ee..e079c77 100644
--- a/src/make_cue_matrix.jl
+++ b/src/make_cue_matrix.jl
@@ -457,6 +457,107 @@ function make_combined_cue_matrix(
     cue_obj_train, cue_obj_val
 end
 
+
+"""
+    make_cue_matrix_from_CFBS(features::Vector{Vector{T}};
+                              pad_val::T = 0.,
+                              ncol::Union{Missing,Int}=missing) where {T}
+
+Create a cue matrix from a vector of feature vectors (usually CFBS vectors).
+It is expected (though of course not necessary) that the vectors have varying lengths. They are consequently padded on the right with the provided `pad_val`.
+
+# Obligatory arguments
+- `features::Vector{Vector{T}}`: vector of vectors containing C-FBS features
+
+# Optional arguments
+- `pad_val::T = 0.`: Value with which the feature vectors will be padded
+- `ncol::Union{Missing,Int}=missing`: Number of columns of the C matrix. If not set, will be set to the maximum number of features
+
+# Examples
+```julia
+C = JudiLing.make_cue_matrix_from_CFBS(features)
+```
+"""
+function make_cue_matrix_from_CFBS(features::Vector{Vector{T}};
+                                    pad_val::T = 0.,
+                                    ncol::Union{Missing,Int}=missing) where {T}
+
+    m = findmax([length(v) for v in features])[1]
+
+    if !ismissing(ncol)
+        if ncol >= m
+            m = ncol
+        else
+            error("ncol is set to be smaller than the maximum number of features")
+        end
+    end
+
+    n = length(features)
+
+    C = fill(pad_val, (n, m))
+    for i in 1:n
+        C[i, 1:length(features[i])] = features[i]
+    end
+
+    return C
+end
+
+
+"""
+    make_combined_cue_matrix_from_CFBS(features_train::Vector{Vector{T}},
+                                       features_test::Vector{Vector{T}};
+                                       pad_val::T = 0.,
+                                       ncol::Union{Missing,Int}=missing) where {T}
+
+
+Create cue matrices from two vectors of feature vectors (usually CFBS vectors).
+It is expected (though of course not necessary) that the vectors have varying lengths. They are consequently padded on the right with the provided `pad_val`.
+The cue matrices are set to have to the size of the maximum number of feature values in `features_train` and `features_test`.
+
+# Obligatory arguments
+- `features_train::Vector{Vector{T}}`: vector of vectors containing C-FBS features
+- `features_test::Vector{Vector{T}}`: vector of vectors containing C-FBS features
+
+# Optional arguments
+- `pad_val::T = 0.`: Value with which the feature vectors will be padded
+- `ncol::Union{Missing,Int}=missing`: Number of columns of the C matrices. If not set, will be set to the maximum number of features in `features_train` and `features_test`
+
+# Examples
+```julia
+C_train, C_test = JudiLing.make_combined_cue_matrix_from_CFBS(features_train, features_test)
+```
+"""
+function make_combined_cue_matrix_from_CFBS(features_train::Vector{Vector{T}},
+                                            features_test::Vector{Vector{T}};
+                                            pad_val::T = 0.,
+                                            ncol::Union{Missing,Int}=missing) where {T}
+
+    m = findmax([length(v) for v in [features_train;features_test]])[1]
+
+    if !ismissing(ncol)
+        if ncol >= m
+            m = ncol
+        else
+            error("ncol is set to be smaller than the maximum number of features")
+        end
+    end
+
+    n_train = length(features_train)
+    n_test = length(features_test)
+
+    C_train = fill(pad_val, (n_train, m))
+    for i in 1:n_train
+        C_train[i, 1:length(features_train[i])] = features_train[i]
+    end
+
+    C_test = fill(pad_val, (n_test, m))
+    for i in 1:n_test
+        C_test[i, 1:length(features_test[i])] = features_test[i]
+    end
+
+    return C_train, C_test
+end
+
 """
     make_ngrams(tokens, grams, keep_sep, sep_token, start_end_token)
 
diff --git a/test/eval_tests.jl b/test/eval_tests.jl
index 2538b3c..934c847 100644
--- a/test/eval_tests.jl
+++ b/test/eval_tests.jl
@@ -171,6 +171,61 @@ end
 
 end
 
+@testset "token_accuracy" begin
+    S = [[1. 2. 3.]
+         [4. 5. 6.]
+         [7. 8. 9.]]
+
+    S2 = [[1. 2. 3.]
+          [0. 0.1 0.1]
+          [9. 8. 7.]]
+
+    @test JudiLing.eval_SC(S, S2) ≈ 0.3333
+    @test JudiLing.eval_SC(S, S2, freq=[1,1,1]) ≈ 0.3333
+    @test JudiLing.eval_SC(S, S2, freq=[1,2,3]) ≈ 0.1667
+    @test JudiLing.eval_SC(S, S2, freq=[3,2,1]) ≈ 0.5
+    @test JudiLing.eval_SC(S, S, freq=[1,2,3]) == JudiLing.eval_SC(S, S)
+
+    S3 = [7. 8. 9.]
+
+    @test JudiLing.eval_SC(S, S2, S3) ≈ 0.3333
+    @test JudiLing.eval_SC(S, S2, S3, freq=[1,1,1]) ≈ 0.3333
+    @test JudiLing.eval_SC(S, S2, S3, freq=[1,2,3]) ≈ 0.1667
+    @test JudiLing.eval_SC(S, S2, S3, freq=[3,2,1]) ≈ 0.5
+    @test JudiLing.eval_SC(S2, S2, S3, freq=[1,2,3]) == JudiLing.eval_SC(S2, S2)
+
+    S4 = [[1. 2. 3.]
+         [6. 5. 4.]
+         [7. 8. 9.]]
+
+    S5 = [[1. 2. 3.]
+          [9. 3. 1.]
+          [0. 0.1 0.1]]
+
+    S6 = [6. 5. 4.]
+
+    @test JudiLing.eval_SC(S4, S5) ≈ 0.6667
+    @test JudiLing.eval_SC(S4, S5, S6) ≈ 0.3333
+    @test JudiLing.eval_SC(S4, S5, freq=[1,1,1]) ≈ 0.6667
+    @test JudiLing.eval_SC(S4, S5, S6, freq=[1,1,1]) ≈ 0.3333
+    @test JudiLing.eval_SC(S4, S5, S6, freq=[1,2,3]) ≈ 0.1667
+    @test JudiLing.eval_SC(S4, S5, S6, freq=[3,2,1]) ≈ 0.5
+
+    data = DataFrame("Word"=>["a", "b", "a"])
+    @test JudiLing.eval_SC(S, S2, data, "Word") ≈ 0.6667
+    @test JudiLing.eval_SC(S, S2, data, "Word", freq=[1,1,1]) ≈ 0.6667
+    @test JudiLing.eval_SC(S, S2, data, "Word", freq=[1,3,2]) ≈ 0.5
+    @test JudiLing.eval_SC(S, S2, data, "Word", freq=[3,1,2]) ≈ round(5/6, digits=4)
+    @test JudiLing.eval_SC(S2, S2, data, "Word", freq=[3,1,2]) ≈ 1.0
+
+    data = DataFrame("Word"=>["a", "b", "c"])
+    data2 = DataFrame("Word"=>["b"])
+    @test JudiLing.eval_SC(S4, S5, S6, data, data2, "Word") ≈ 0.6667
+    @test JudiLing.eval_SC(S4, S5, S6, data, data2, "Word", freq=[1,1,1]) ≈ 0.6667
+    @test JudiLing.eval_SC(S4, S5, S6, data, data2, "Word", freq=[3,2,1]) ≈ round(5/6, digits=4)
+    @test JudiLing.eval_SC(S4, S5, S6, data, data2, "Word", freq=[1,2,3]) ≈ 0.5
+end
+
 @testset "accuracy_comprehension" begin
     latin = DataFrame(
         Word = ["ABC", "BCD", "CDE", "BCD"],
diff --git a/test/make_cue_matrix_tests.jl b/test/make_cue_matrix_tests.jl
index ff60021..9415e0c 100644
--- a/test/make_cue_matrix_tests.jl
+++ b/test/make_cue_matrix_tests.jl
@@ -123,3 +123,61 @@ end
         keep_sep = false,
     )
 end
+
+
+@testset "CFBS" begin
+
+    features = [[0., 0.1, 1.0, 3.2, 3.1], [-1., 0.5], [-0.]]
+
+    C = JudiLing.make_cue_matrix_from_CFBS(features)
+    target_C = [[0. 0.1 1.0 3.2 3.1]
+                [-1. 0.5 0. 0. 0.]
+                [-0. 0. 0. 0. 0.]]
+    @test C == target_C
+
+    C = JudiLing.make_cue_matrix_from_CFBS(features, ncol = 6)
+    target_C = [[0. 0.1 1.0 3.2 3.1 0.]
+                [-1. 0.5 0. 0. 0. 0.]
+                [-0. 0. 0. 0. 0. 0.]]
+    @test C == target_C
+
+    C = JudiLing.make_cue_matrix_from_CFBS(features, pad_val = 1.)
+    target_C = [[0. 0.1 1.0 3.2 3.1]
+                [-1. 0.5 1. 1. 1.]
+                [-0. 1. 1. 1. 1.]]
+    @test C == target_C
+
+    @test_throws TypeError JudiLing.make_cue_matrix_from_CFBS(features, pad_val = 1)
+
+    @test_throws ErrorException JudiLing.make_cue_matrix_from_CFBS(features, ncol=4)
+
+    features2 = [[0., 0.1, 1.0, 3.2, 3.1, 4.05]]
+
+    C1, C2 = JudiLing.make_combined_cue_matrix_from_CFBS(features, features2)
+    target_C1 = [[0. 0.1 1.0 3.2 3.1 0.]
+                [-1. 0.5 0. 0. 0. 0.]
+                [-0. 0. 0. 0. 0. 0.]]
+    target_C2 = [0. 0.1 1.0 3.2 3.1 4.05]
+    @test C1 == target_C1
+    @test C2 == target_C2
+
+    C1, C2 = JudiLing.make_combined_cue_matrix_from_CFBS(features, features2, ncol = 7)
+    target_C1 = [[0. 0.1 1.0 3.2 3.1 0. 0.]
+                [-1. 0.5 0. 0. 0. 0. 0.]
+                [-0. 0. 0. 0. 0. 0. 0.]]
+    target_C2 = [0. 0.1 1.0 3.2 3.1 4.05 0.]
+    @test C1 == target_C1
+    @test C2 == target_C2
+
+    C1, C2 = JudiLing.make_combined_cue_matrix_from_CFBS(features, features2, pad_val = 1.)
+    target_C1 = [[0. 0.1 1.0 3.2 3.1 1.]
+                [-1. 0.5 1. 1. 1. 1.]
+                [-0. 1. 1. 1. 1. 1.]]
+    target_C2 = [0. 0.1 1.0 3.2 3.1 4.05]
+    @test C1 == target_C1
+    @test C2 == target_C2
+
+    @test_throws TypeError JudiLing.make_combined_cue_matrix_from_CFBS(features, features2, pad_val = 1)
+
+    @test_throws ErrorException JudiLing.make_combined_cue_matrix_from_CFBS(features, features2, ncol=4)
+end