diff --git a/Examples/CS1_1.jl b/Examples/CS1_1.jl
index ad6363e..4bd07f4 100644
--- a/Examples/CS1_1.jl
+++ b/Examples/CS1_1.jl
@@ -4,7 +4,7 @@ using JLD2
 
 ## CS1_1  - Case-study example: Monitoring flight actuator faults
 #           No measurements of control surface angles are used
-println("Case study CS1_1 with Fig8.2 and Fig8.3")
+println("Case study CS1_1 with Fig8.2(fig1) and Fig8.3(fig2)")
 
 ## Part 1 - Model setup
 # load matrices of the aircraft multiple-model SYSACM, 
@@ -132,11 +132,11 @@ display([NormRu NormRd NormRfmRfnom])
 # generate figures
 y, tout, x = stepresp(Rftilde,10);
 include("Fig8_2.jl")
-Fig8_2 = f
+Fig8_2 = fig1
 
 y = [stepresp(Rtot[i][:,indf],10)[1] for i in 1:N];
 include("Fig8_3.jl")
-Fig8_3 = f
+Fig8_3 = fig2
 
 end  # module
 
diff --git a/Examples/CS1_2.jl b/Examples/CS1_2.jl
index cf40e90..f9d67c4 100644
--- a/Examples/CS1_2.jl
+++ b/Examples/CS1_2.jl
@@ -121,7 +121,7 @@ display([NormRu NormRd NormRfmRfnom])
 y = [stepresp(Rtot[i][:,indf],2)[1] for i in 1:N];
 tout = Vector(0:0.02:2.)
 include("Fig8_4.jl")
-Fig8_4 = f
+Fig8_4 = fig
 
 end  # module
 using Main.CS1_2
diff --git a/Examples/CS2_1.jl b/Examples/CS2_1.jl
index b63eb26..60e383e 100644
--- a/Examples/CS2_1.jl
+++ b/Examples/CS2_1.jl
@@ -5,7 +5,7 @@ using FaultDetectionTools, DescriptorSystems, LinearAlgebra,
 
 ## CS2_1  - Case-study example: Monitoring air data sensor faults
 #           Robust least order LTI synthesis
-println("Case study CS2_1 with Fig8.5, Fig8.6 and Fig8.7")
+println("Case study CS2_1 with Fig8.5(fig1), Fig8.6(fig2) and Fig8.7(fig3)")
 
 ## Part 1 - Model setup
 # load matrices of the aircraft multiple-model SYSACSM and 
@@ -81,7 +81,7 @@ end
 
 tout = Vector(0:0.1:10.)
 include("Fig8_5.jl")
-Fig8_5 = f
+Fig8_5 = fig1
 
 ##  Part 5 - Assesment of synthesis results for nominal synthesis
 
@@ -99,7 +99,7 @@ end
 
 tout = Vector(0:0.1:10.)
 include("Fig8_6.jl")
-Fig8_6 = f
+Fig8_6 = fig2
 
 
 ##  Part 6 - Assesment of synthesis results for mean value gains
@@ -196,7 +196,7 @@ display(norms)
 
 tout = Vector(0:0.1:10.)
 include("Fig8_7.jl")
-Fig8_7 = f
+Fig8_7 = fig3
 
 end  # module
 using Main.CS2_1
diff --git a/Examples/CS2_2.jl b/Examples/CS2_2.jl
index abc2673..e277ecb 100644
--- a/Examples/CS2_2.jl
+++ b/Examples/CS2_2.jl
@@ -86,7 +86,7 @@ display(Nintp)
 
 tout = Vector(0:0.1:10.)
 include("Fig8_8.jl")
-Fig8_8 = f
+Fig8_8 = fig
 
 end  # module
 using Main.CS2_2
diff --git a/Examples/Ex5_16c.jl b/Examples/Ex5_16c.jl
index 2a13550..ae9b36b 100644
--- a/Examples/Ex5_16c.jl
+++ b/Examples/Ex5_16c.jl
@@ -54,7 +54,7 @@ yl = y .- yerr; yu = y .+ yerr  # lower and upper bounds
 
 # plot step responses
 include("Fig5_2.jl")
-Fig5_2 = f
+Fig5_2 = fig
 
 end
 using Main.Ex5_16c
diff --git a/Examples/Ex6_1c.jl b/Examples/Ex6_1c.jl
index 09d5fc5..dd2ae2f 100644
--- a/Examples/Ex6_1c.jl
+++ b/Examples/Ex6_1c.jl
@@ -2,7 +2,7 @@ module Ex6_1c
 using FaultDetectionTools, DescriptorSystems, LinearAlgebra, Test
 
 # Example 6.1c - Solution of an EMDP
-println("Example 6.1c with Fig6.1 and Fig6.2")
+println("Example 6.1c with Fig6.1(fig1) and Fig6.2(fig2)")
 
 # Lateral aircraft model without faults
 A = [-.4492 0.046 .0053 -.9926;
diff --git a/Examples/Ex7_1.jl b/Examples/Ex7_1.jl
index f47ef23..e3bc532 100644
--- a/Examples/Ex7_1.jl
+++ b/Examples/Ex7_1.jl
@@ -5,7 +5,7 @@ using DescriptorSystems, Polynomials
 Makie.inline!(false)
 
 # Example 7.1 - Illustrating polynomial root sensitivity 
-println("Example 7.1 with Fig7_1")
+println("Example 7.1 with Fig7.1")
 
 
 Makie.inline!(false)
diff --git a/Examples/Ex7_1a.jl b/Examples/Ex7_1a.jl
index 2b667e0..ba96007 100644
--- a/Examples/Ex7_1a.jl
+++ b/Examples/Ex7_1a.jl
@@ -6,7 +6,7 @@ Makie.inline!(false)
 
 
 # Example 7.1 - Illustrating high precision polynomial root computation 
-println("Example 7.1 with alternative Fig7_1a")
+println("Example 7.1 with alternative Fig7.1a")
 
 
 pexact = BigFloat.(Vector(-25.:1.:-1.))
diff --git a/Examples/Fig5_2.jl b/Examples/Fig5_2.jl
index 55b895a..f5b5f38 100644
--- a/Examples/Fig5_2.jl
+++ b/Examples/Fig5_2.jl
@@ -4,15 +4,16 @@
 #    yl   - the lower bounds of step responses
 #    yu   - the upper bounds of step responses 
 #    tout - the time samples  
-using CairoMakie, LaTeXStrings
+using Makie, CairoMakie, LaTeXStrings
+Makie.inline!(false)
 
 title_u = [ latexstring("From: \$f_1\$")  latexstring("From: \$f_2\$")  latexstring("From: \$u_1\$")  latexstring("From: \$u_2\$")]
 ylabel_r = [ latexstring("To: \$r_1\$")  latexstring("To: \$r_2\$") ]
 
 ns, p, m = size(y)
-f = Figure(resolution = (800, 600))
+fig = Figure(resolution = (800, 600))
 
-axs = [Axis(f[row, col]) for row in 1:p, col in 1:m]
+axs = [Axis(fig[row, col]) for row in 1:p, col in 1:m]
 
 for row in 1:p
     for col in 1:m
@@ -31,17 +32,17 @@ for row in 1:p
     end
 end
 
-Label(f[0, :], text = "Step Responses",  fontsize = 20, 
+Label(fig[0, :], text = "Step Responses",  fontsize = 20, 
                font = "TeX Gyre Heros Bold", valign = :bottom, 
                padding = (0, 0, -10, 0))
-Label(f[end+1, :], text = "Time (seconds)", font = "TeX Gyre Heros Bold",
+Label(fig[end+1, :], text = "Time (seconds)", font = "TeX Gyre Heros Bold",
                    valign = :top,padding = (0, 0, 5, -10))
-Label(f[0:end, 0], text = "Residuals", font = "TeX Gyre Heros Bold", rotation = pi/2, 
+Label(fig[0:end, 0], text = "Residuals", font = "TeX Gyre Heros Bold", rotation = pi/2, 
                    valign = :center, padding = (0, -20, 0, 0))
 axs[1,1].yticks = [0,0.5,1,1.5]
 axs[2,1].yticks = [0,0.5,1,1.5]
 
-f
+fig
 
 # comment out next line to save plot
-# save("Fig5_2.pdf", f, resolution = (800, 600))
+# save("Fig5_2.pdf", fig, resolution = (800, 600))
diff --git a/Examples/Fig8_2.jl b/Examples/Fig8_2.jl
index 48c053c..6c50610 100644
--- a/Examples/Fig8_2.jl
+++ b/Examples/Fig8_2.jl
@@ -11,9 +11,9 @@ yticks = [[-20,0,20],[-20,0,20],[-20,0,20],[0,10,20],[0,1,2],[0,1,2]]
 yhighs = [20,20,20,20,2,2]
 ylows = [-20,-20,-20,-1,-0.1,-0.1]
 ns, p, m = size(y)
-f = Figure(resolution = (800, 600))
+fig1 = Figure(resolution = (800, 600))
 
-axs = [Axis(f[row, col]) for row in 1:p, col in 1:m]
+axs = [Axis(fig1[row, col]) for row in 1:p, col in 1:m]
 
 for row in 1:p
     for col in 1:m
@@ -33,10 +33,10 @@ for row in 1:p
     end
 end
 
-Label(f[end+1, :], text = "Time (seconds)", font = "TeX Gyre Heros Bold",
+Label(fig1[end+1, :], text = "Time (seconds)", font = "TeX Gyre Heros Bold",
                    valign = :top,padding = (0, 0, 5, -10))
 
-f
+fig1
 
 # comment out next line to save plot
-#save("Fig8_2.pdf", f, resolution = (800, 600))
+#save("Fig8_2.pdf", fig1, resolution = (800, 600))
diff --git a/Examples/Fig8_3.jl b/Examples/Fig8_3.jl
index 392891c..d0fe328 100644
--- a/Examples/Fig8_3.jl
+++ b/Examples/Fig8_3.jl
@@ -13,9 +13,9 @@ yticks = [[-50,0,50],[-50,0,50],[-50,0,50],[-50,0,50],[-10,0,10],["   0","   1",
 yhighs = [50,50,50,51,10,2]
 ylows = [-50,-51,-51,-50,-10,-0.2]
 ns, p, m = size(y[1])
-f = Figure(resolution = (800, 600))
+fig2 = Figure(resolution = (800, 600))
 
-axs = [Axis(f[row, col]) for row in 1:p, col in 1:m]
+axs = [Axis(fig2[row, col]) for row in 1:p, col in 1:m]
 
 for row in 1:p
     for col in 1:m
@@ -37,10 +37,10 @@ for row in 1:p
     end
 end
 
-Label(f[end+1, :], text = "Time (seconds)", font = "TeX Gyre Heros Bold",
+Label(fig2[end+1, :], text = "Time (seconds)", font = "TeX Gyre Heros Bold",
                    valign = :top,padding = (0, 0, 5, -10))
 
-f
+fig2
 
 # comment out next line to save plot
-#save("Fig8_3.pdf", f, resolution = (800, 600))
+#save("Fig8_3.pdf", fig2, resolution = (800, 600))
diff --git a/Examples/Fig8_4.jl b/Examples/Fig8_4.jl
index 444bdf6..c5821a6 100644
--- a/Examples/Fig8_4.jl
+++ b/Examples/Fig8_4.jl
@@ -9,9 +9,9 @@ N = length(y)
 title_u = reshape([ latexstring("From: \$f_$i\$")  for i in 1:mf],1,mf)
 ylabel_r = reshape([ latexstring("To: \$r_$i\$")  for i in 1:mf],1,mf)
 ns, pp, mm = size(y[1])
-f = Figure(resolution = (800, 600))
+fig = Figure(resolution = (800, 600))
 
-axs = [Axis(f[row, col]) for row in 1:pp, col in 1:mm]
+axs = [Axis(fig[row, col]) for row in 1:pp, col in 1:mm]
 
 for row in 1:pp
     for col in 1:mm
@@ -34,10 +34,10 @@ for row in 1:pp
     end
 end
 
-Label(f[end+1, :], text = "Time (seconds)", font = "TeX Gyre Heros Bold",
+Label(fig[end+1, :], text = "Time (seconds)", font = "TeX Gyre Heros Bold",
                    valign = :top,padding = (0, 0, 5, -10))
 
-f
+fig
 
 # comment out next line to save plot
-# save("Fig8_4.pdf", f, resolution = (800, 600))
+# save("Fig8_4.pdf", fig, resolution = (800, 600))
diff --git a/Examples/Fig8_5.jl b/Examples/Fig8_5.jl
index bbc312d..5b64bbb 100644
--- a/Examples/Fig8_5.jl
+++ b/Examples/Fig8_5.jl
@@ -9,9 +9,9 @@ N = length(y)
 title_u = [ L"From: $f_1$"  L"From: $f_2$"  L"From: $u_1$"  L"From: $u_2$" L"From: $u_3$" L"From: $d_1$"  L"From: $d_2$"]
 ylabel_r = [ latexstring("To: \$r_1\$")  latexstring("To: \$r_2\$") ]
 ns, pp, mm = size(y[1])
-f = Figure(resolution = (800, 500))
+fig1 = Figure(resolution = (800, 500))
 
-axs = [Axis(f[row, col]) for row in 1:pp, col in 1:mm]
+axs = [Axis(fig1[row, col]) for row in 1:pp, col in 1:mm]
 inputs = [mu+md .+ (1:mf); 1:mu+md]
 for row in 1:pp
     for col in 1:mm
@@ -34,10 +34,10 @@ for row in 1:pp
     end
 end
 
-Label(f[end+1, :], text = "Time (seconds)", font = "TeX Gyre Heros Bold",
+Label(fig1[end+1, :], text = "Time (seconds)", font = "TeX Gyre Heros Bold",
                    valign = :top,padding = (0, 0, 5, -10))
 
-f
+fig1
 
 # comment out next line to save plot
-# save("Fig8_5.pdf", f, resolution = (800, 500))
+# save("Fig8_5.pdf", fig1, resolution = (800, 500))
diff --git a/Examples/Fig8_6.jl b/Examples/Fig8_6.jl
index aeddf43..1c6923e 100644
--- a/Examples/Fig8_6.jl
+++ b/Examples/Fig8_6.jl
@@ -12,9 +12,9 @@ yticks = [["   -1", "0", "1", "2"],[-50,0,50,100]]
 yhighs = [2,100]
 ylows = [-1,-50]
 ns, pp, mm = size(y[1])
-f = Figure(resolution = (800, 500))
+fig2 = Figure(resolution = (800, 500))
 
-axs = [Axis(f[row, col]) for row in 1:pp, col in 1:mm]
+axs = [Axis(fig2[row, col]) for row in 1:pp, col in 1:mm]
 inputs = [mu+md .+ (1:mf); 1:mu+md]
 for row in 1:pp
     for col in 1:mm
@@ -37,10 +37,10 @@ for row in 1:pp
     end
 end
 
-Label(f[end+1, :], text = "Time (seconds)", font = "TeX Gyre Heros Bold",
+Label(fig2[end+1, :], text = "Time (seconds)", font = "TeX Gyre Heros Bold",
                    valign = :top,padding = (0, 0, 5, -10))
 
-f
+fig2
 
 # comment out next line to save plot
-# save("Fig8_6.pdf", f, resolution = (800, 500))
+# save("Fig8_6.pdf", fig2, resolution = (800, 500))
diff --git a/Examples/Fig8_7.jl b/Examples/Fig8_7.jl
index 90ee1a9..5b16e60 100644
--- a/Examples/Fig8_7.jl
+++ b/Examples/Fig8_7.jl
@@ -9,9 +9,9 @@ N = length(y)
 title_u = [ L"From: $f_1$"  L"From: $f_2$"  L"From: $u_1$"  L"From: $u_2$" L"From: $u_3$" L"From: $d_1$"  L"From: $d_2$"]
 ylabel_r = [ latexstring("To: \$r_1\$")  latexstring("To: \$r_2\$") ]
 ns, pp, mm = size(y[1])
-f = Figure(resolution = (800, 500))
+fig3 = Figure(resolution = (800, 500))
 
-axs = [Axis(f[row, col]) for row in 1:pp, col in 1:mm]
+axs = [Axis(fig3[row, col]) for row in 1:pp, col in 1:mm]
 inputs = [mu+md .+ (1:mf); 1:mu+md]
 for row in 1:pp
     for col in 1:mm
@@ -34,10 +34,10 @@ for row in 1:pp
     end
 end
 
-Label(f[end+1, :], text = "Time (seconds)", font = "TeX Gyre Heros Bold",
+Label(fig3[end+1, :], text = "Time (seconds)", font = "TeX Gyre Heros Bold",
                    valign = :top,padding = (0, 0, 5, -10))
 
-f
+fig3
 
 # comment out next line to save plot
-# save("Fig8_7.pdf", f, resolution = (800, 500))
+# save("Fig8_7.pdf", fig3, resolution = (800, 500))
diff --git a/Examples/Fig8_8.jl b/Examples/Fig8_8.jl
index 5b6750f..ca26d0b 100644
--- a/Examples/Fig8_8.jl
+++ b/Examples/Fig8_8.jl
@@ -9,9 +9,9 @@ N = length(y)
 title_u = [ L"From: $f_1$"  L"From: $f_2$"  L"From: $u_1$"  L"From: $u_2$" L"From: $u_3$" L"From: $d_1$"  L"From: $d_2$"]
 ylabel_r = [ latexstring("To: \$r_1\$")  latexstring("To: \$r_2\$") ]
 ns, pp, mm = size(y[1])
-f = Figure(resolution = (800, 500))
+fig = Figure(resolution = (800, 500))
 
-axs = [Axis(f[row, col]) for row in 1:pp, col in 1:mm]
+axs = [Axis(fig[row, col]) for row in 1:pp, col in 1:mm]
 inputs = [mu+md .+ (1:mf); 1:mu+md]
 for row in 1:pp
     for col in 1:mm
@@ -34,10 +34,10 @@ for row in 1:pp
     end
 end
 
-Label(f[end+1, :], text = "Time (seconds)", font = "TeX Gyre Heros Bold",
+Label(fig[end+1, :], text = "Time (seconds)", font = "TeX Gyre Heros Bold",
                    valign = :top,padding = (0, 0, 5, -10))
 
-f
+fig
 
 # comment out next line to save plot
-# save("Fig8_8.pdf", f, resolution = (800, 500))
+# save("Fig8_8.pdf", fig, resolution = (800, 500))
diff --git a/Project.toml b/Project.toml
index 9b3a56b..e8734ed 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "FaultDetectionTools"
 uuid = "fbfc5fa0-0005-11ec-19bb-b3a492551a0b"
 authors = ["Andreas Varga <varga.andreas@gmail.com>"]
-version = "1.1.1"
+version = "1.2.0"
 
 [deps]
 Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"
diff --git a/README.md b/README.md
index cb8a240..284c8af 100644
--- a/README.md
+++ b/README.md
@@ -51,7 +51,7 @@ using FaultDetectionTools
 cd(joinpath(pkgdir(FaultDetectionTools), "Examples"))
 include("CS2_1.jl")
 ````
-_Note:_ For the execution of the test examples and case study examples, the packages **Measurements**, **Makie**, **CairoMakie**, **LaTeXStrings**, **JLD2** and **Optim** are also required and must be additionally installed. 
+_Note:_ For the execution of the test examples and case study examples, the packages **Measurements**, **GenericLinearAlgebra**, **Makie**, **CairoMakie**, **LaTeXStrings**, **JLD2** and **Optim** are also required and must be additionally installed. 
 
 ## About
 
diff --git a/ReleaseNotes.md b/ReleaseNotes.md
index aa3733a..e1a49d5 100644
--- a/ReleaseNotes.md
+++ b/ReleaseNotes.md
@@ -1,5 +1,9 @@
 # Release Notes
 
+## Version 1.2.0
+
+New analysis function `mdgenspec` and changed functionality and improvements in `mdsspec` and `mdspec`.  
+
 ## Version 1.1.1
 
 Version bump to use Makie 0.19. Bugs fixed in MDObjects.jl.
diff --git a/docs/src/MDanalysis.md b/docs/src/MDanalysis.md
index 9d508c1..c5aeae7 100644
--- a/docs/src/MDanalysis.md
+++ b/docs/src/MDanalysis.md
@@ -1,9 +1,11 @@
 # Analysis of model detection synthesis models
 
+* **[`mdgenspec`](@ref)**  Generation of achievable model detection specifications.
 * **[`mddist`](@ref)**  Computation of distances between component models.
 * **[`mddist2c`](@ref)**  Computation of pairwise distances between two sets of component models.
 
 ```@docs
+mdgenspec
 mddist
 mddist2c
 ```
diff --git a/docs/src/index.md b/docs/src/index.md
index 7aba29a..fd89ee0 100644
--- a/docs/src/index.md
+++ b/docs/src/index.md
@@ -71,6 +71,7 @@ The available functions in the current version of the `FaultDetectionTools.jl` p
 
 **Analysis of model detection synthesis models**
 
+* **[`mdgenspec`](@ref)**  Generation of achievable model detection specifications.
 * **[`mddist`](@ref)**  Computation of distances between component models.
 * **[`mddist2c`](@ref)**  Computation of pairwise distances between two sets of component models.
 
diff --git a/src/FDIanalysis.jl b/src/FDIanalysis.jl
index 02f35b6..dbfd4e2 100644
--- a/src/FDIanalysis.jl
+++ b/src/FDIanalysis.jl
@@ -194,7 +194,7 @@ The keyword argument `sdeg = β` specifies a prescribed stability degree `β` fo
 generated candidate filters, such that the real parts of
 filters poles must be less than or equal to `β`, in the continuous-time case, and 
 the magnitudes of filter poles must be less than or
-equal to `β`, in the discrete-time case. If `sdeg = missing` then no then no stabilization is performed if and `FDFreq = missing`.
+equal to `β`, in the discrete-time case. If `sdeg = missing` then no stabilization is performed if `FDFreq = missing`.
 If `sdeg = missing` and `FDFreq = freq`, then the following default values are employed : `β = -0.05`, in continuous-time case, and  `β = 0.95`, 
 in discrete-time case. 
 
diff --git a/src/FaultDetectionTools.jl b/src/FaultDetectionTools.jl
index d554c73..1098c27 100644
--- a/src/FaultDetectionTools.jl
+++ b/src/FaultDetectionTools.jl
@@ -22,7 +22,7 @@ export efdsyn, efdisyn, efdbasesel, afdsyn, afdredsyn, afdisyn, afdbasesel, emms
 export fdhinfminus, fdhinfmax, binmat2dec, dec2binmat, fditspec_, fdisspec_, fdiscond_
 
 export MDMModel, MDModel, mdmodset, MDFilter, MDFilterIF, mdIFeval
-export mddist, mddist2c
+export mdgenspec, mddist, mddist2c
 export mdspec, mdsspec, mdperf, mdmatch, mdgap
 export emdsyn, amdsyn, emdbasesel
 
diff --git a/src/MDanalysis.jl b/src/MDanalysis.jl
index 73b57cf..2da4e22 100644
--- a/src/MDanalysis.jl
+++ b/src/MDanalysis.jl
@@ -1,3 +1,159 @@
+"""
+    S = mdgenspec(sysm::MDMModel; emtest = false, fast = true, MDtol, MDGainTol, MDfreq, 
+                  sdeg, atol, atol1, atol2, atol3, rtol) 
+                           
+Generate for a given multiple synthesis model `sysm::MDMModel` the achievable binary structure matrix `S` 
+of the internal form corresponding to a set of stable model detection filters, chosen such that the `i`-th residual is 
+insensitive to the `i`-th model for all control and disturbance inputs. 
+If `N` is the number of component models, then `S` is an `N×N` binary matrix, whose `(i,i)`-th 
+element is zero for `i = 1, ..., N`, 
+and its `(i,j)`-th element is nonzero if there exists a model detection filter such that the corresponding 
+`i`-th residual is sensitive to the `j`-th model for certain control inputs (if `emtest = false`) or 
+for certain control and disturbance inputs (if `emtest = true`).   
+   
+`MDFreq = freq` specifies a vector of real frequency values or a scalar real frquency value
+for strong model detectability checks (default: `MDFreq = missing`).
+
+`MDtol = tol1` specifies the threshold `tol1` for model detectability checks
+   (default: `tol1 = 0.0001`).
+
+`MDGainTol = tol2` specifies the threshold `tol2` for strong model detectability checks
+   (default: `tol2 = 0.01`). 
+
+The keyword argument `sdeg = β` specifies a prescribed stability degree `β` for the poles of the internally 
+generated candidate filters, such that the real parts of
+filters poles must be less than or equal to `β`, in the continuous-time case, and 
+the magnitudes of filter poles must be less than or
+equal to `β`, in the discrete-time case. If `sdeg = missing` then no stabilization is performed if `FDFreq = missing`.
+If `sdeg = missing` and `FDFreq = freq`, then the following default values are employed : `β = -0.05`, in continuous-time case, and  `β = 0.95`, 
+in discrete-time case. 
+
+The rank determinations in the performed reductions
+are based on rank revealing QR-decompositions with column pivoting 
+if `fast = true` (default) or the more reliable SVD-decompositions if `fast = false`.
+
+The keyword arguments `atol1`, `atol2` and `rtol`, specify, respectively, 
+the absolute tolerance for the nonzero elements of the state-space matrices (see below)
+`Ai`, `Bi`, `Ci` and `Di`, the absolute tolerance for the nonzero elements of `Ei`,   
+and the relative tolerance for the nonzero elements of all above matrices.  
+The default relative tolerance is `ni*ϵ`, where `ϵ` is the working machine epsilon 
+and `ni` is the orders of the system `sysm.sys[i]`. 
+The keyword argument `atol` can be used 
+to simultaneously set `atol1 = atol` and `atol2 = atol`. 
+
+   
+_Method:_ For a multiple synthesis model `sysm` containing `N` stable component models, 
+the `i`-th component system `sysm.sys[i]` has a descriptor realization of the form 
+`sysm.sys[i] = (Ai-λEi,Bi,Ci,Di)` and the corresponding input-output form is
+
+      yi = Gui(λ)*u + Gdi(λ)*di + Gwi(λ)*wi + Gvi(λ)*vi 
+
+with the Laplace- or Z-transformed plant outputs `yi`, control inputs `u`, 
+disturbance inputs `di`, noise inputs `wi`, and auxiliary inputs `vi`,  
+and with `Gui(λ)`, `Gdi(λ)`, `Gwi(λ)` and 
+`Gvi(λ)`, the corresponding transfer function matrices. 
+
+To generate the achievable structure matrix `S`, it is assumed that `N` model detection filters are used, 
+`Q_i(λ)`, for `i = 1, ..., N`, where the `i`-th model detection filter `Q_i(λ)` generates the i-th residual
+
+            r_i = Q_i(λ)*| y |
+                         | u |               
+
+and is determined such that it fulfills the nullspace synthesis condition
+ 
+            Q_i(λ)*|Gui(λ) Gdi(λ)| = 0.
+                   |  I     0    |
+
+To generate the elements of the `i`-th row of `S`, the corresponding internal forms of the `i`-th filter 
+are determined for `j = 1, ..., N` as                    
+
+       | Ruij(λ) Rdij(λ) | := Q_i(λ)*| Guj(λ) Gdj(λ) | 
+                                     |   I      0    |
+                   
+and `S[i,j]` is set as follows: 
+   - if `MDfreq = missing` and `emtest = false`, 
+     then `S[i,j] = 1` if `Ruij(λ) ̸= 0` or `S[i,j] = 0` if `Ruij(λ) = 0`;
+   - if `MDfreq = missing` and `emtest = true`, 
+     then `S[i,j] = 1` if `[Ruij(λ) Rdij(λ)] ̸= 0` or `S[i,j] = 0` if `[Ruij(λ) Rdij(λ)] = 0`;
+   - if `MDfreq = freq` and `emtest = false`, 
+     then `S[i,j] = 1` if `Ruij(λs) ̸= 0` for any `λs ∈ freq` or 
+          `S[i,j] = 0` if `Ruij(λs) = 0` for all `λs ∈ freq` ;
+   - if `MDfreq = freq` and `emtest = true`, 
+     then `S[i,j] = 1` if `[Ruij(λs) Rdij(λs)] ̸= 0` for any `λs ∈ freq` or 
+          `S[i,j] = 0` if `[Ruij(λs) Rdij(λs)] = 0` for all `λs ∈ freq` .
+
+
+_References:_
+
+[1] A. Varga, Solving Fault Diagnosis Problems - Linear Synthesis Techniques. 
+              Springer Verlag, 2017; sec. 5.2.
+
+[2] A. Varga, Least order fault and model detection using multi-models. 
+       IEEE CDC'09, Shanghai, China, 2009.
+"""
+function mdgenspec(sysm::MDMModel; sdeg::Union{Real,Missing} = missing, emdtest::Bool = false,
+                      MDtol::Real = 0.0001, MDGainTol::Real = 0.01, MDfreq::Union{AbstractVector{<:Real},Real,Missing} = missing, 
+                      atol::Real = zero(Float64), atol1::Real = atol, atol2::Real = atol, 
+                      rtol::Real = (size(sysm.sys[1].A,1)+1)*eps(Float64)*iszero(max(atol1,atol2)), 
+                      fast::Bool = true) 
+   N = length(sysm.sys)   
+   T1 = Float64
+   # check stability of input model
+   isstable(sysm) || error("the input multiple model must be stable")
+
+   disc = (sysm.sys[1].Ts != 0);  # system type (continuous- or discrete-time)
+  
+   # decode options
+    
+   strongMD = !ismissing(MDfreq)
+   strongMD && (freq = isa(MDfreq,Vector) ? MDfreq : [MDfreq]) 
+  
+   # set default stability degree
+   sdegdefault = disc ? 0.95 : -0.05
+   
+   ismissing(sdeg) && (sdeg = sdegdefault)  # set desired stability degree to default value
+  
+ 
+   # decode input information
+   mu = sysm.mu  
+   md = sysm.md 
+   m = mu .+ md      # total number of inputs
+   sdegNS = strongMD ? sdegdefault : missing   
+   S = falses(N,N)  
+   for i = 1:N
+       # solve the exact MD problem for the i-th system (with noise inputs)
+       # 
+       # compute a minimal left nullspace basis Q1i for the i-th system
+      # form [ Gu Gd Gf Gw Gaux; I 0 0 0 0] 
+      syse = [sysm.sys[i][:,1:m[i]]; eye(T1,mu,m[i])];
+      #
+      # compute a left nullspace basis Q = Q1 of G1 = [Gu Gd; I 0] = 0 and
+      # obtain QR = [ Q R ], where R = [ Rw Ra] = Q*[Gw Ga;0 0 0]
+      qtemp = glnull(syse; atol1, atol2, rtol, fast, sdeg = sdegNS)[1] 
+      # check solvability conditions
+      size(qtemp,1) == 0 && (@warn "empty nullspace basis for the $i-th model")
+
+
+      #strongMD && (qtemp = glcf(qtemp; atol1, atol2, rtol, fast, sdeg = sdegNS)[1])
+      for j = 1:N
+         # check j-th model detectability
+         if i != j
+            mj = emdtest ? m[j] : mu
+            temp = gir(qtemp*[sysm.sys[j][:,1:mj]; eye(T1,mu,mj)]; atol1, atol2, rtol)
+            if strongMD
+               S[i,j] = any(fdisspec_(temp, freq; FDGainTol = MDGainTol, atol1, atol2, rtol = 0, fast, stabilize = false)[1])
+            else
+               S[i,j] = any(fditspec_(temp; atol1, atol2, rtol, FDtol = MDtol))
+            end
+         end
+      end
+   end
+
+   return S
+   
+   # end MDGENSPEC
+end
+
 """
      mddist2c(sysm1, sysm2; MDfreq, cdinp = false, distance = "nugap", 
               rtolinf = 0.00001, offset, atol, atol1 = atol, atol2 = atol, rtol = 0, fast = true) -> (dist, fpeak)
diff --git a/src/MDperformance.jl b/src/MDperformance.jl
index 1586852..570a851 100644
--- a/src/MDperformance.jl
+++ b/src/MDperformance.jl
@@ -32,12 +32,9 @@ function mdspec(sysR::MDFilterIF{T};  cdinp::Bool = false,
     mu = sysR.mu
     cdinp && (md = sysR.md)
     for j = 1:N
+        inpj = cdinp ? (1:mu+md[j]) : (1:mu)
         for i = 1:M
-            if cdinp 
-               S[i,j] = !iszero(sysR.sys[i,j][:,1:mu+md[j]]; atol1, atol2, rtol)
-            else
-               S[i,j] = !iszero(sysR.sys[i,j][:,1:mu]; atol1, atol2, rtol)
-            end
+            S[i,j] = !iszero(sysR.sys[i,j][:,inpj]; atol1, atol2, rtol)
         end         
     end
     return S
@@ -67,8 +64,8 @@ Let the `(i,j)`-th component filter `sysR.sys[i,j]` have the input-output form
 with the Laplace- or Z-transformed residual output `rij`, control inputs `u`, 
 disturbance inputs `dj`, noise inputs `wj`, and auxiliary inputs `vj`,  
 and with `Ruij(λ)`, `Rdij(λ)`, `Rwij(λ)` and `Rvij(λ)`, the corresponding transfer function matrices. 
-Then, `S[i,j] = 1` if `Ruij(λs)` is nonzero for all `λs ∈ Ω` and `cdinp = false` (default), or
-if `[Ruij(λs) Rdij(λs)]` is nonzero for all `λs ∈ Ω` and `cdinp = true`. Otherwise, `S[i,j] = 0`.
+Then, `S[i,j] = 1` if `Ruij(λs)` is nonzero for any `λs ∈ Ω` and `cdinp = false` (default), or
+if `[Ruij(λs) Rdij(λs)]` is nonzero for any `λs ∈ Ω` and `cdinp = true`. Otherwise, `S[i,j] = 0`.
 
 `MDGainTol = tol` specifies an absolute  threshold `tol` for the nonzero magnitudes of 
 the frequency response gains (default: `tol = 0.01`). 
@@ -90,19 +87,14 @@ function mdsspec(sysR::MDFilterIF{T}, freq::Union{AbstractVector{<:Real},Real} =
                   MDGainTol::Real = 0.01, atol::Real = zero(float(real(T))), atol1::Real = atol, atol2::Real = atol, 
                   rtol::Real = 0, fast::Bool = true) where T
     isa(freq,Vector) ? tfreq = freq : tfreq = [freq] 
-    lfreq = length(tfreq);
     M, N = size(sysR.sys)
     S = trues(M,N)
     mu = sysR.mu
     cdinp && (md = sysR.md)
     for j = 1:N
+        inpj = cdinp ? (1:mu+md[j]) : (1:mu)
         for i = 1:M
-            if cdinp 
-               Sc = fdisspec_(sysR.sys[i,j][:,1:mu+md[j]], tfreq; FDGainTol = MDGainTol, block = true, atol1, atol2, rtol, fast)[1]
-            else
-               Sc = fdisspec_(sysR.sys[i,j][:,1:mu], tfreq; FDGainTol = MDGainTol, block = true, atol1, atol2, rtol, fast)[1]
-            end
-            S[i,j] = lfreq == 1 ? any(Sc) : all([any(view(Sc,:,:,k)) for k = 1:lfreq])  
+            S[i,j] = any(fdisspec_(sysR.sys[i,j][:,inpj], tfreq; FDGainTol = MDGainTol, block = true, atol1, atol2, rtol, fast)[1])
         end         
     end
     return S
diff --git a/src/types/MDObjects.jl b/src/types/MDObjects.jl
index b772851..b8c156e 100644
--- a/src/types/MDObjects.jl
+++ b/src/types/MDObjects.jl
@@ -491,3 +491,15 @@ gpole(Q::MDFilter{T}; kwargs...) where T = gpole.(Q.sys;  kwargs...)
 isstable(Q::MDFilter{T}; kwargs...) where T = all(isstable.(Q.sys;  kwargs...))
 isstable(sysm::MDMModel; kwargs...) = all(isstable.(sysm.sys;  kwargs...))
 
+
+function Base.getproperty(sys::MDFilter, d::Symbol)  
+    if d === :outputs
+        return 1:sys.ny
+    elseif d === :controls
+        return sys.ny+1:sys.ny+sys.mu
+    else
+        getfield(sys, d)
+    end
+end
+Base.propertynames(sys::MDFilter) =
+    (fieldnames(typeof(sys))..., :outputs, :controls)
diff --git a/test/test_ammsyn.jl b/test/test_ammsyn.jl
index 341d7da..558a204 100644
--- a/test/test_ammsyn.jl
+++ b/test/test_ammsyn.jl
@@ -262,7 +262,7 @@ Q, R, info = afdisyn(sysf, SFDI; smarg =-3, sdeg = -3)
 
 Mr = FDFilterIF([R.sys[1][:,R.faults]; R.sys[2][:,R.faults]];mf);
 
-@time QM, RM, info1 = ammsyn(sysf, Mr; atol = 1.e-7, reltol = 1.e-5, sdeg = -3);
+@time QM, RM, info1 = ammsyn(sysf, Mr; atol = 1.e-7, reltol = 1.e-4, sdeg = -3);
 
 R1 = fdIFeval(QM,sysf); # form Q*[Gu Gd Gf;I 0 0];
 
diff --git a/test/test_emdsyn.jl b/test/test_emdsyn.jl
index 9b2630c..f468445 100644
--- a/test/test_emdsyn.jl
+++ b/test/test_emdsyn.jl
@@ -40,6 +40,12 @@ display(sysm)
 @time distgap3c, fpeak3c = mddist2c(sysm, sysm, distance = "Inf")
 @test norm(distgap3-distgap3c,Inf) < 1.e-7 
 
+@time S = mdgenspec(sysm)
+@test S == ([0  1  1  1; 1  0  1  1; 1  1  0  1; 1  1  1  0] .> 0)
+@time S0 = mdgenspec(sysm,MDfreq = 0)
+@test all(S0 .== false)
+@test S == mdgenspec(sysm,MDfreq = [0,1])
+
 @time Q, R, info = emdsyn(sysm; sdeg = -15, poles = [-20]); info.MDperf
 display(Q)
 display(R)