Find discontinuities and lift up or pull down the part after the discontinuity

Hi everyone i have a singal with discontinuities.Does anyone know how can i find the discontinuities and at those points lift up or pull down the remaining signal?

[Jan]

Please post some real data. The diagram is not clear enough. What is the width of the jumps? Is this one time step? What about using a low pass filter?

[lena kappa]

Sorry!

This is one of my datasets (my x axes is just 1:800)

193

196

197

200

203

205

207

209

211

214

215

217

218

219

221

221

222

222

222

222

221

220

219

218

216

213

212

209

205

203

199

197

192

188

185

180

175

173

168

165

161

156

154

150

146

144

141

140

137

135

134

133

133

133

133

134

135

136

138

140

143

145

148

150

154

157

160

164

167

170

173

175

178

181

182

184

185

186

186

186

185

184

183

180

177

175

171

168

163

158

154

148

141

137

130

126

119

112

108

102

95

91

86

82

77

73

70

67

64

62

61

60

59

59

60

61

62

63

66

67

70

74

76

80

84

86

91

93

97

101

103

107

111

113

116

118

120

123

124

126

128

128

130

130

131

132

132

132

132

132

132

131

131

130

130

129

128

127

127

126

125

124

123

122

122

121

120

120

119

119

119

118

119

119

119

120

120

121

122

123

125

126

128

130

132

134

136

139

142

143

147

150

151

154

156

158

161

162

165

167

167

169

169

170

171

171

171

170

170

168

167

166

164

162

160

157

155

151

149

145

142

139

135

131

129

125

122

119

115

113

109

106

104

101

100

98

96

94

93

92

92

91

91

91

91

92

92

94

94

96

97

99

100

102

103

105

107

108

110

111

113

114

115

116

116

117

117

117

117

117

116

116

115

113

112

111

109

107

105

102

101

98

96

93

90

88

85

81

79

76

74

71

68

66

63

60

58

56

53

52

50

49

47

46

45

44

43

43

43

43

43

44

45

47

48

50

52

54

57

60

62

66

70

72

76

79

83

88

90

95

99

102

106

108

112

116

118

121

124

126

128

130

132

134

135

136

138

138

139

139

140

140

140

140

140

140

140

140

139

139

138

138

137

136

135

134

133

132

131

130

128

127

125

124

122

120

119

117

115

114

113

111

110

109

108

107

107

106

106

105

105

106

106

106

107

107

108

108

109

109

110

110

111

111

111

110

110

110

109

108

108

106

105

104

103

101

100

98

97

95

94

93

92

91

91

90

91

91

92

92

93

95

96

99

102

104

107

109

113

117

119

123

127

129

133

137

139

142

144

147

150

151

153

155

155

156

156

157

157

156

156

154

153

152

151

149

146

145

142

140

138

135

133

130

128

126

124

121

120

118

116

115

113

112

112

111

110

110

109

110

110

110

110

111

111

112

112

113

114

114

115

115

115

116

115

116

115

115

114

114

113

112

110

109

107

106

104

101

100

97

95

93

91

89

87

85

83

81

80

78

77

76

75

75

74

74

74

75

75

75

77

78

78

80

81

82

84

85

86

88

89

90

91

93

94

95

96

97

97

98

99

99

100

100

101

101

102

102

102

103

103

103

104

104

105

105

106

107

107

108

109

110

110

111

113

113

115

115

117

118

119

120

121

122

123

124

125

126

126

127

128

128

128

128

128

128

128

128

127

126

125

124

123

121

120

118

116

114

112

110

107

105

103

100

98

96

94

93

91

89

88

88

87

87

87

87

88

89

90

92

94

95

98

100

103

107

109

113

117

120

124

126

130

134

137

140

144

146

149

151

154

156

157

159

160

161

162

163

163

163

163

163

162

161

160

160

158

157

156

154

152

150

148

147

144

142

140

138

136

134

132

130

128

126

124

122

120

118

116

115

113

112

111

109

108

107

106

106

105

105

104

104

104

104

104

104

105

105

105

106

107

107

108

109

110

111

111

112

112

113

113

114

114

114

114

113

113

112

111

110

110

108

107

106

104

103

101

99

98

95

94

92

90

89

87

85

84

82

81

80

79

78

77

77

77

77

77

77

78

79

80

82

83

85

87

90

93

95

99

103

105

109

112

117

121

124

129

133

136

140

143

147

150

152

155

157

159

160

161

162

which gives me this

[Jan]

I do noit see the steps in your diagram.

Answer 1

[Jan]

x = linspace(1, 10*pi, 200);
y = sin(x);
y(60:90) = y(60:90) + 1;
y(120:180) = y(120:180) - 1;
yOrig = y;

limit = 0.2;
dy    = [0, diff(y)];
jump  = strfind(abs(dy) > limit, [false, true, false]);
for ijump = 1:numel(jump)
   k = jump(ijump);
   y(k+1:end) = y(k+1:end) - dy(k+1);
end
figure
axes('NextPlot', 'add')
plot(x, yOrig, '-r', x, y, 'bo')
y2 = lowpass(y, 30, 200);
plot(x, y2, 'c+')

[lena kappa]

Could you explain to me how this code works?

[Jan]

@lena kappa: The code searchs for steps greater than the limit. Then it deletes the height of the step from the rest of the vector.

lowpass is a low pass filter: A sharp step is a hight frequency signal. Lettting only low frequencies pass, removes jumps.

[lena kappa]

Hi @Jan i modified your code this way

x_0 = [1 800];
y_0 = [1 501.06];

y = improfile(img101,x_0,y_0)
yOrig = y;


limit = 0.2;
dy    = [0, diff(y)];
jump  = strfind(abs(dy) > limit, [false, true, false]);
for ijump = 1:numel(jump)
   k = jump(ijump);
   y(k+1:end) = y(k+1:end) - dy(k+1);
end
figure
axes('NextPlot', 'add')
plot(x, yOrig, '-r', x, y, 'bo')
y2 = lowpass(y, 30, 200);
plot(x, y2, 'c+')

and i get the next error:

Error using horzcat
Dimensions of arrays being concatenated are not consistent.
Error in lowpass_discontinuities_simple (line 34)
dy    = [0, diff(y)];

do you know why?

[lena kappa]

all my code is :

img101 = imread('c1=100.800.1.0.1.png');

img_rot= imrotate(img101,32.06);

   prof_img_rot(1:256,i)=img_rot(212+i,400:655)';


x_0 = [1 800];
y_0 = [1 501.06];

y = improfile(img101,x_0,y_0)
yOrig = y;


limit = 0.2;
dy    = [0, diff(y)];
jump  = strfind(abs(dy) > limit, [false, true, false]);
for ijump = 1:numel(jump)
   k = jump(ijump);
   y(k+1:end) = y(k+1:end) - dy(k+1);
end
figure
axes('NextPlot', 'add')
plot(x, yOrig, '-r', x, y, 'bo')
y2 = lowpass(y, 30, 200);
plot(x, y2, 'c+')

[Jan]

The command

dy    = [0, diff(y)];

is working, if y is a row vector. For a column vector replace the comma by a semicolon for a horizontal concatenation:

dy    = [0; diff(y)];

[lena kappa]

Yes yes that's it! Thanks @Jan

[lena kappa]

@Jan how did you choose the numbers 30 and 200 in

y2 = lowpass(y, 30, 200);

[Jan]

@lena kappa: See the documentation of lowpass:

y = lowpass(x,fpass,fs) specifies that x has been sampled at a rate of fs hertz. fpass is the passband frequency of the filter in hertz.

The result is stable over a wide range of parameters.

Answer 2

[Image Analyst]

See a similar discussion started/asked by Matt J:

https://www.mathworks.com/matlabcentral/answers/895202-fit-a-transform-to-remove-piecewise-discontinuities?s_tid=srchtitle

[lena kappa]

Hello and thanks for your response hiw do you choose a and b in your code?

% Create uncorrupted data as a noisy sine wave.
x = linspace(pi/4,pi/2,60);
Y = sin(x) + randn(size(x))/80;
% plot(X, Y, 'b--s', 'LineWidth', 2);
% hold on
% Lower a chunk of it.
a=1.1; %unknown
b=-0.5; %unknown
I = 1 <= x & x <= 1.3; %KNOWN

[Image Analyst]

@lena kappa That's not my code. You'd have to ask the author of the code.

[lena kappa]

But you uploaded this code in the link you sent above.

clc;    % Clear the command window.
close all;  % Close all figures (except those of imtool.)
clear;  % Erase all existing variables. Or clearvars if you want.
workspace;  % Make sure the workspace panel is showing.
format long g;
format compact;
fontSize = 15;
% Create uncorrupted data as a noisy sine wave.
x = linspace(pi/4,pi/2,60);
Y = sin(x) + randn(size(x))/80;
% plot(X, Y, 'b--s', 'LineWidth', 2);
% hold on
% Lower a chunk of it.
a=1.1; %unknown
b=-0.5; %unknown
I = 1 <= x & x <= 1.3; %KNOWN
Yc = Y;
Yc(I) = a*Y(I)+b; %corrupted signal
% Plot original signal.
subplot(2, 1, 1);
plot(x, Yc, '--s', 'LineWidth', 2);
grid on;
xlabel('x', 'FontSize', fontSize);
ylabel('y', 'FontSize', fontSize);
title('Original Signal', 'FontSize', fontSize);
% Now we have data and we can begin.
%======================================================
% Get a smoothed signal
windowSize = 31;
ym = movmax(Yc, windowSize);
% Plot it.
hold on;
plot(x, ym, '-', 'LineWidth', 2);
% Find the difference between smoothed and actual
diffSignal = abs(ym - Yc);
plot(x, diffSignal, '-', 'LineWidth', 2);
% The lowered signal is where the difference is more than 3
% Get the good data.
threshold = 0.3;
yline(threshold, 'Color', 'g', 'LineWidth', 3);
goodIndexes = diffSignal < threshold;
badIndexes = diffSignal >= threshold;
xGood = x(goodIndexes);
yGood = Yc(goodIndexes);
% Interpolate over the bad section.
ym = interp1(xGood, yGood, x);
% Plot the interpolated section
plot(x, ym, 'm-');
% Update the difference signal.
diffSignal = abs(ym - Yc);
% Get average difference over ths bad section
amountToLift = mean(diffSignal(badIndexes))
% List bad data.
yRepaired = Yc; % Initialize
yRepaired(badIndexes) = Yc(badIndexes) + amountToLift;
% Plot repaired signal.
subplot(2, 1, 2);
plot(x, yRepaired, 'g.-', 'LineWidth', 2, 'MarkerSize', 20);
grid on;
xlabel('x', 'FontSize', fontSize);
ylabel('y', 'FontSize', fontSize);
title('Repaired Signal', 'FontSize', fontSize);
g = gcf;
g.WindowState = 'maximized';

[lena kappa]

i think i understood it! sorry to bother you again

[Image Analyst]

I guess those were just some parameters to cause a chunk of the waveform to shift by a reasonable amount.