texts:
I am interested in how should a listing for a program look like that will generate a sine of f=50hz, using the next arguments for any moment in time.

float phis=162*PI/180;//initial phase of signal
////////////////////////////////////////////////////////////////////////
//  generates a sinusoid of dwLength in dwData
//      and maintains phis value updated
////////////////////////////////////////////////////////////////////////
void rc(short* dwDataOut, DWORD dwLength){
float omegaf=50/(float)dwLength/PI*90/100;
DWORD dwi;
//use phis
    for (dwi=0;dwi<dwLength;dwi++)
        dwDataOut[dwi]=4450*sin(dwi*omegaf+phis)+200;

//update wave phase
    phis-=PI*(1-50/(float)dwLength)+PI/2;
    while(phis<0)
        phis+=2*PI;
}

The program I came up works correctly with a small marj of error only for f=50hz, because right from the start it appears with a small error in getting phis coefficient.

By applying this procedure in a program I'll make use of a cleaner signal from a sound wave, and obtain a signal that will be overlapped as antiphase signal to a more complex 50hz wave, for longer period of time, and so eliminating this harmonic. I’m interested in a more mathematical method to find the value for phis, and probably the arguments of sine function.

(is this function resulting in a true output for f=50hz? what are the errors that can be adjusted?)


texts: How do I calculate the noise variance in this problem?

texts:
In a paper about radar, the transmitted signal expressed in low-pass equivalent form is given by $\sqrt{E}s(t)$, where $E$ is the transmitted energy, $s(t)$ satisfies
$$
\int_T | s(t) |^2 dt = 1,
$$
and $T$ is the pulse duration.

Here
$$
s(t) = \sqrt{\frac{1}{T}} \exp(j2\pi f_c t)
$$
and the received signal is
$$
r(t) = \alpha s(t) + n(t)
$$
and the SNR in this paper is
$$
\frac{E |\alpha|^2}{\sigma^2}
$$
where $\alpha$ is the reflectivity, and $\sigma^2$ is the variance of the zero-mean Gaussian noise process $n(t)$.

So if I want to generate $r(t)$ in Matlab with SNR = 10, I sample $s(t)$ with a sample period $T_S$ so I get
$$
s(nT_s) = \sqrt{\frac{1}{T}} \exp(j2\pi f_c nT_S)
$$

My question is:  how do I calculate the noise variance, $\sigma^2$ from the above so that I can calculate:
$$
r(nT_S) = \alpha s(nT_s) + \sigma {\tt randn(size(}s{\tt ))}
$$

??


texts:
Q1:  Consider a discrete-time band-pass filter (BPF) that uses only two poles as

$$ H(z) = G \frac{1}{\big(1-r \, e^{j \omega_0} z^{-1}\big)\big(1-r \, e^{-j \omega_0} z^{-1}\big)} $$

where this filter has a center frequency of $\omega_0$.  If this BPF was designed to have a very narrow bandwidth by selecting $r$ to be very close to unity, show that the 3dB bandwidth of the BPF filter is approximately given by

$$ BW_{3dB}=2(1-r) $$

(Hint:  Use a geometrical approach on the unit circle.)


texts: What is the between spatial frequency of an image and pixel size of the sensor?

texts:
In my lecture notes about the sampling of an image I've written that:

Since pixels have a finite dimension the spatial frequency response is attenuated before the "ideal" Nyquist frequency
In order to have a better freq. response, the sensor resolution should increase
Also the sensor dimension should increase otherwise the pixel simply gets smaller and the problem remains

I really can't understand the last sentence. Why should a small pixel size be negative for the spatial frequency response? I'm probably misunderstanding (or I've misunderstood while listening to the lecture) the whole thing because in my mind it should be the opposite.


texts: Unexpected results for fourier resampling (Python 3)

texts:
I am trying to resample a signal using the fourier method in python using np.fft.rfft() and np.fft.irfft(). (I'm not using scipy.signal.resample() because in my final application I need to do more to my signal in the frequency domain than what can be accomplished by the window argument of scipy.signal.resample().) I have reduced the problems I see to this simple test case (ipython notebook):

%pylab inline

N = 16
P = 8

S = [0.49, 1.68, 0.78, 0.05, 0.21, 0.29, 1.34, 1.17, 16.73,
     48.78, 16.90, 0.62, 0.40, 1.60, 0.86, 0.57]
assert len(S) == N

figure(figsize=(14, 5))
plot(np.arange(N), S, '.')

dft = np.fft.rfft(S)
plot(np.arange(N), np.fft.irfft(dft), '-')

dft *= P
dft.resize(P*N/2+1)
plot(np.arange(P*N)/P, np.fft.irfft(dft), '-')

xlim(6.5, 11.5)
ylim(-5, 55);

And here is what I get back:

As you can see, the result looks almost correct, but is nevertheless a bit off. (The resampled red curve does not go through the original sample points, even though there are sample points at the same positions in time as in the original signal.)

The greeen line shows that everything is still ok before extending and scaling the dft array. So the problem seems to be in the actual resampling:

dft *= P
dft.resize(P*N/2+1)
plot(np.arange(P*N)/P, np.fft.irfft(dft), '-')

I don't see how this result is mathematically possible because other than scaling of time and amplitude axis, adding zeroes to the high-frequency end of the dft array should only add zero terms to the fourier series expansion of the signal. So the signal should still pass through the original sample values.

But the small error I see does not seem to be related to scaling: sometimes the resampled waveform lags behind the original (t=8) and sometimes it is faster (t=10), sometimes it is larger in magnitude (t=7,9,11) and sometimes it is smaller (t=8,10).

I'm playing with this for hours now and can't find what is wrong with my code.. Maybe I just don't see it because I'm already staring at it for so long now..

(see np.fft documentation for the definitions used by the python FFT implementations)

Edit (in response to comment): When I replace S with an impulse at t=9, or a fast impulse train, I get the following:


texts: Get an input signal from LTI system output, find impulse response

texts:
I just study signal processing and I have some questions.

How to find impulse response $h[n]$ sequence of the system if I have
$x[n]$ and $y[n]$ (input and output) sequences? Hope there is a formula
to do that.
How to find $x[n]$ input sequence if I have $h[n]$ and $y[n]$?
Is it possible to find $x[n]$ if I just have set of $y[n]$ values, but I
don't know impulse response $h[n]$? How to do that if so?

Maybe you can also suggest something for further reading about this topics.


texts: Understand about registration between spatial prior and original image

texts:
I am working in level set/active contour method using spatial prior information . The main papers that I am focusing:"Automatic segmentation of neonatal images using convex optimization and coupled level sets"-2011. It is very nice papers.

However, I am confusing about some points in the paper. In the paper, the author said that "For a given subject image, by using HAMMER algorithm to warp the atlas image into the subject image space, we can derive spatial prior information priori for WM, GM and CSF, respectively, for the subject image".As my understand, we have a original image (need to segment), a template image, and a spatial prior image (range [0,1]). Using HAMMER method, you will transform spatial prior image in to original image. Then, original image is fixed and prior image is moving. Is it right? Or vice versa, the spatial prior image is fixed and original image is moving, which is correct? Because the main reason of warping step is that we want to find a corresponding prior value of a given input pixel in the spatial prior map.

When I read some tool such as SPM, I saw that original image often fixed and spatial prior information is moving. Is it correct?

Thanks in advance


texts: Power Spectral Density Via Rank 1 Update

texts:
I know that I can find the autocorrelation matrix of a series of finite length sequences via rank-1 updates using the relation:

$\mathbf{R}[k] = \frac{1}{k}\mathbf{R}[k-1] + \frac{1}{k+1}\mathbf{r}[k]\mathbf{r}^{H}[k]$

Where $k$ is the discrete time index and $\mathbf{r}$ is the autocorrelation of the current signal sequence. The true autocorrelation can then be found via the first column of the matrix (see https://en.wikipedia.org/wiki/Autocorrelation_matrix)

So I imagine I can do the same for the PSD, correct? Namely:

$\mathbf{S}[k] = \frac{1}{k}\mathbf{S}[k-1] + \frac{1}{k+1}\mathbf{s}[k]\mathbf{s}^{H}[k]$

Where $\mathbf{s}$ is the fourier spectrum of the current frame. The PSD will then be found as the first column of $\mathbf{S}$. Correct?

Is there a way I can transform $\mu: \mathbf{R} \rightarrow \mathbf{S}$ directly, that is, without having to compute the PSD of every sequence frame?


texts:
I would like to understand a concept of homodyne architecture of I-Q demodulator and how the wave looks like at each step.

It is a bandpass Signal:

Step 1: multiplication with cos and sin wave.  I represent cos and sin wave like Dirac function. For cos is even, sin is odd.

Step 2: filter with LPF

Step 3: summation

Will be it correct?


texts: How does the quantization error generate noise?

texts: I'm learning about sampling and DSP on my own. I have a hard time to understand how the quantization error results in noise. I think I miss a fundamental understanding but can't tell what it is. So how does the quantization error generate noise?


texts:
I'm doing an exercise where, given a data sample taken from a physical measure, I have to estimate the parameters of the measure where samples $x_i$ are described by the equation:
$$
x_i = A \sin(2\pi f t_i) + Y_i + Z_i
$$
where:

$A \sin(2\pi f t)$ is a bias introduced by the measure
$Y ∼ \operatorname{Logis}(\mu, s)$ is the R.V. describing the population
$Z ∼ N(0, \sigma_Z)$ is Gaussian noise

I have to estimate the five parameters $A,f,\mu, s, \sigma_Z$.

I've found $\mu$ which simply corresponds to the sample mean $\bar{X}$.
To estimate $f$ I think I can analyze the autocorrelation but I'm not sure on how to do it.

My main problem is on finding the variances of $Y$ and $Z$. I'm trying to find them using a system with two equations, however I only have:
$$s^2 = \sigma_Y^2 + \sigma_Z^2$$
which indicates that the estimated variance $s^2$, since $Y$ and $Z$ are independent R.V, is just the sum of the two variances.

Is there any other relation between those two variances I don't spot?


texts:
Say input-output of a system is defined as:

$$
x[n] \longrightarrow x[nM]
$$

then what will be the output of $x[n-1]$?

will it be:

\begin{align}
x[n-1] \longrightarrow  &x\left[(n-1)M\right] = x\left[nM - M\right]\\
&\textrm{or}\\
x[n-1] \longrightarrow  &x\left[nM - 1\right]
\end{align}

And if someone says shift a signal by 1, if the signal is $x[2n-1]$ then the result will be

$$x[2n-1-1]\quad \textrm{or}\quad x[2(n-1)-1]\quad ?$$

I'm really confused between these. Can anyone make me understand it, please?


texts: Find the period of a signal with the DTFT plot

texts:
I have an exercise and I'm struggling to resolve it. Here it is :

My problem is about the DTFT. I've always been taught that we use DTFT for infinite-lenght signal that are not periodic (if the signal were periodic I would have used DFS)... So, how am I supposed to find the period of x[n] and y[n] ? I must have missed or misunderstood something. Can somebody explain it to me ? (obviously without giving the explicit answer). Thanks !


texts: why butterworth bandpass filter changes signal in the passband?

texts:
I want to apply a butterworth bandpass filter to EEG signals and it's not working so I tried with a simple signal:

t=[0:0.01:10];
s3= 9*sin(2*pi*4*t);

Still there is a problem(probably obvious) that i cannot detect. :/
Here is how I design and apply the filter with passband [3 5]Hz:

%Filtering
Fs=100;
Wp1=3;
Wp2=5;
[n,Wn]=buttord([Wp1 Wp2]/(Fs/2), [Wp1-0.5 Wp2+0.5]/(Fs/2), 4, 40);
[z, p, k] = butter(n, Wn);
[sos,g] = zp2sos(z,p,k);
Xfilt_s3= filtfilt(sos, g, s3);

This is the initial signal and the signal after butterworth application in time and frequency domain.

I expected that the signal wouldn't change after the filtering, because its frequency is in the passband.
Also, why is this amplitude change occurring?
Lastly, since I want to apply it at EEG signals, I would like to know if filtfilt function is appropriate for EEG filtering.
Thank you! :)


texts:
in my attempts to synthesize a plucked electric guitar string (using supercollider's digital waveguide, which is based on the Karplus-Strong algorithm) I came to a dead end and decided to look that the Cepstrums of a sampled string and my synthesized string, to nail down the audible differences. They look like this (A-string - 110Hz)

My questions are:

The sampled signal has something in the right half, which the synthesized signal does not. I cannot see how to interpret this as "excitation". So what does it stand for?
Is it safe to say, that the spectrum of the synthesized spectrum is more "harmonic" than the one of the sampled signal? When I look at the spectra themselves it does appear this way. The sampled spectrum (left) looks more "noisy")

Could it be that the regularity (and the dull sound) of the synthesized signal is a consequence of the digital waverguide itself. The way I understand Karplus-Strong, after some time only the fundamental and its harmonics survive, which would indeed lead to a "harmonic" spectrum.
If that is the case, what could be the cause of the irregularties of a real plucked string. IOW: what does a real string have, which Karplus-Strong doesn't have?
As a side note: what does the sharp peak at the end of the Cepstrum stand for? Is it caused by the 1 second window?


texts:
I am trying to estimate the average amplitude of some signal with frequency 6 Hz, sampled at ~300 Hz. See figures for a part of the signal and its dft calculated using matlab.
I estimate the average amplitude of the signal by calculating its root mean square,
>> rms(signal)

ans =

    5.1651

However, when looking at the signal (see first figure), it is very clear the mean amplitude is slightly greater than ~5, and is probably closer to ~8. Summing the peak amplitude in the frequency domain and its first few aliases results in the same value (~5).
What am I doing wrong here? How to properly estimate the signal average amplitude from its dft?


texts: Difference between DC component and zero frequency component of signal

texts:
We know that Fourier Transform of a signal exists if it is absolutely integrable
and it exists for periodic signals if impulse functions are allowed.

If we consider the fourier transform of $\text{rect}(t)$ , we get $\text{sinc}(f)$ in frequency domain. Which has got zero frequency component as $1$. But we all know, DC value of $\text{rect}(t)$ is zero.

My question is:

If a signal has got a zero frequency component  in frequency domain ,There must be DC value in time domain. But why there is no DC value in case of $\text{rect}(t)$ in time domain?
What is the difference between DC component and zero frequency component?


texts: Reconstruct a signal via continuous wavelet transform (CWT)

texts:
I want to use the continuous wavelet transform (CWT) to modify the amplitude of a signal in the frequency domain (in a specific frequency band) and then reconstruct the signal. I want to do it by scaling the wavelet coefficients in that band. Is that the right way?

Thanks in advance for your help.


texts: FIR filter design with time and frequency domain constraints

texts:
I deal with the problem of the synthesis of the filter with the requirements of both the frequency and
time characteristics.

Now there is a filter that meets the requirements.
Roughly speaking the requirements are as follows:

Frequency response

FIR Low Pass

Fpass ~ 0.0(3) (* pi rad/sample)
Fstop ~ 0.1(6) (* pi rad/sample)
bandwidth irregularity of ~ 1e-3 (maybe ~ 1e-2) dB
~ -40 dB suppression

Time requirements

Overshoot in step response ~ 4% (transition time ~ 1e-2 ms)

I would like to know about the possibilities of such a synthesis (I’m with filterDesigner in the matlab failed, tried free tools - Iowa Hills)
or at least about the transformation of the finished filter into a narrower
(Fpass = 0.025, Fstop = 0.125) (Matlab can do this, but only IIR is obtained).


texts:
I need to visualize some data as gray-scale images after contrast enhancement. The processing algorithm must be built on few samples, and then it will be applied on unseen images (suppose, the [min, max] range is given).

I plan to:

Use sample images to estimate actual histogram H1.
Define desired histogram H2.
Build a transformation T (look-up table) that transforms H1 into H2.
Apply T on unseen images.

Question: How to build H2? Are there studies in regards to the optimal gray-scale histogram for Human eye?

Below are two examples of different transformations, applied on the same image:


texts: Given a Horizontal Derivative Kernel What Would Be the Output Image?

texts: Hi, im solving this problem and everything went well until i went and solved the position p(5,3) and i got the  answer of 0. However, the solution guide  lists the solution  as 150 which is pointed at in blue. I don't understand this since the pixel at the location is 0, meaning we don't have to calculate the solution however the solution says otherwise. Thanks for your time.


texts:
I've been working for a while on 3D building reconstruction from old photos (c. 1904).  In all of these cases, the cameras and buildings are long gone, making standard camera calibration impossible.  To even get an affine reconstruction, the only thing I have to rely on (that I know of) is vanishing points.
Even ignoring image measurement error, vanishing points have a few problems:

I'm at the mercy of the craftspeople who built the buildings.  In many cases, the allegedly parallel lines don't converge very well to a single point.  And the actual orthogonality of the vanishing points is unknown.
I'm at the mercy of the photographer.  In most cases, the photos are nearly 2-point perspectives.  There is usually enough tilt that I can estimate the third vanishing point, but it's a poor estimate, and it's often a couple of orders of magnitude further away than the other two, which I assume further degrades the estimate.

Initially, I had hoped that I could get a rough reconstruction and improve it by manually "tweaking" the vanishing points.  Unfortunately, there's no clear relationship between the reconstruction error and the vanishing point error.  I've looked some at bundle adjustment, but in most cases I only have a few photos that include any particular feature.
Does anyone have any suggestions of how to improve or supplement the vanishing point accuracy?  Here are a couple of things I've tried:

Calculate the trifocal tensor from three images, and use it to check the accuracy of a vanishing point by transferring the same vanishing point from the other two images.
Manually correct the reconstruction, find the (3D) homography between it and the original reconstruction, and use this to correct the cameras.

Do these seem reasonable?  Does anyone have any other/better suggestions?


texts: How to add a particular noise on a signal?

texts:
I have a simple ADC model (an input signal == sampled signal ). The adc should be modelled as with following values:

1.5 V (peak to peak).. meaning --> +1.5 V --> 0 dBFS ... how to do this ?
The noise on the signal should be -120 dBFS/Hz in bandwidth of 2 MHz.. .how to add this exact noise ?
The NFFT= 2^11 and fs= 200 MHz.
The input will be 10 MHz single tone, and then the FFT will be perfomed.

fs = 200e6;

% Adding Noise
n_en=1; % Noise On/Off (1/0)
BW= 2e6;
noise= -120; %dBFS/Hz
SNR=abs(noise + 10*log10(BW));

N= 2^11;
Df_fft = fs/N;
t = (0:N-1)/fs; t = t(:);

f=10e6;

f = Df_fft*round(f/Df_fft);

y = cos(2*pi*f*t);  % I think i need some scaling factor here for  +1.5 V --> 0 dBFS

if n_en==1
    y_noise = awgn(y,SNR,'measured');
else
    y_noise= y;
end

S = 2*fft(y_noise,N)/N;
S = 20*log10(abs(S));
f = (fs/1e6)*linspace(0,1,N);
plot(f,S)

So my questions are what NFFT points to take.?
as that would also determine the noise floor.. ..
what exactly would the noise floor be. ?. how exactly i can see the -120 dBFS/Hz in the graph or what does it means in terms of the graph...

is this the correct way to add noise ?

Thanks

All help would be appreciated.


texts:
I went through the answer Kalman filter in practice and it seems we must know all the first and second order properties of random variables to apply the Kalman filter.

But when I only have a set of data $y_1,\dots,y_n$ and a given model $G$ and $F$ with assumed noise processes, could someone please explain to me how to apply Kalman filter?

Thanks.


texts: Least-squares digital IIR filter design (with arbitrary responses)

texts:
I'm studying the IIR filter design that is described in the book: Algorithms for the constrained design of digital filters with arbitrary phase and magnitude responses.

You can get the code at page 171 (at least the main function), and here is an example of a filter design :

M=4;
N=4;
tau=5;
om=pi*[linspace(0,0.2,20),linspace(0.4,1,60)];
D=[exp(-1i*om(1:20)*tau),zeros(1,60)];
W=ones(1,80);
[b,a,e]=mpiir_l2(M,N,om,D,W,0.98);

This is the design of a pass-band linear phase lowpass filter

There is something I don't understand in this code: exp(-1i*om(1:20)*tau) is a way to create regularly spaced points with constant phase change and magnitude = 1. But I don't understand the parameter tau. I tried to change that parameter and I had totally different results: like a notch instead of a lowpass.

How do you choose that number ?
I have to design my filters according to the user desire but how to guess that parameter?
And also how to you choose the phase angles shift ?

Jeff


texts: FastICA - Purpose of centering and whitening the Data

texts:
I'm learning about ICA, but dont know whats the purpose of centering the data (making each component have zero mean)and whitening the data (not sure what is this for, is it the same as PCA?), I'm following FastICA algorithm for wikipedia
FastICA


texts: How to detect whether the object has triangular shape?

texts:
The task is to detect arrow from medical images. After segmenting, arrow-like objects can be connected with other connected components (noise). In this context, detecting a part of the object that characterizes an arrowhead would be a solution since template based shape descriptor does not perform well.
Please help!

enter image description here


texts: Calculate perimeter and area of Head CT from DICOM

texts: I’m trying to calculate the perimeter and area of the outer lining of the skull from head computed tomography scans. I have these scans as a DICOM, and single slices can be exported as image files. Radiant DICOM viewer can allow you to measure linear distance but not a complex perimeter like the skull. I assume this involves exporting the slice and performing edge analysis in ImageJ, not sure.


texts: Unscented Kalman Filter - Multiple Consecutive Measurement Updates

texts:
In trying to implement an Unscented Kalman Filter (UKF), I have come across the issue of what to do when my measurement signals come in at a different rate than my control inputs, which I use in the prediction step of the UKF. In my implementation, I've tried working around this issue by forming a prediction each time I receive a new measurement, but I've been told that it's also possible and perhaps favorable to only form a prediction when I get a new control input and form a new measurement update when I get a new measurement.

The "standard" Kalman filter algorithm in most textbooks/papers that I've seen do not delve deeply into this issue of different data rates, and they usually just show a simple predict-update loop that occurs at each time step.

I see how it is possible to form multiple consecutive predictions. However, I don't fully understand how one could do an update without a prediction right before it at the same time step, since the computation of the Kalman gain relies on the sigma points passed through the state transition function. Any help would be greatly appreciated!

For reference, I am trying to implement my UKF using the FilterPy Python library: https://filterpy.readthedocs.io/en/latest/_modules/filterpy/kalman/UKF.html#UnscentedKalmanFilter


texts:
How can one correctly downsample a Fourier 2D matrix of complex numbers ? More precisely, if I have a 2D DFT of an image:

x=phantom
X=fft2(x)

How should I do to in the frequency domain to obtain:

XHalf=fft2(imresize(x,0.5))

(And imshow(ifft2(XHalf)) shows the downsampled image of course).

Of course, doing imresize() in the frequency domain does not work.

If you do an imresize(X,0.5), after doing an ifft2(imresize(X,0.5)) the result is nonesense.

I tried simply cropping the Fourier "image": it does give something recognizable but it is clearly not simply that.


texts: What Does It Mean When Cross Correlation Between Two Images Occur at Negative Coordinates

texts:
I am using Matlab normxcorr2 to calculate the cross correlation between several images, something like what has been done here: in MATLAB documentation. There are some of the examples that I get negative values for xoffSet and yoffSet, but I am not sure what it means when these values are negative.
Another issue, is practical since I need to give these xoffSet and yoffSet values as an input to coordinates of rectangle in an image, but since these are negative, the function throws an error. I googled it but no success in finding something that can explain such thing. I was thinking of just putting zero instead of negative values, but I'm not sure if that's a right approach.
I would appreciate if someone can explain about the meaning of such negative values and how I can substitute such negative values with appropriate quantities.


texts:
I'm learning some basic signal processing to use in my research. So I'm trying to solve this combined signal:

$x(t)= \text{sin}(60\pi t+\theta)+  0.5 \text{sin}(120\pi t+\theta)+  0.3 \text{sin}(240\pi t+\theta)$

that are uniformly distributed in range $[0, 2\pi)$. Duration of the combined signal is 20 seconds.

Select a method (e.g. periodogram estimation, the Bartlett method, the Welch method) for estimation of the Power Spectral Density of the combined digital signal And make a parameter selection for the selected method, that provides the following frequency resolution a) 0.05 Hz b) 0.2 Hz c) 0.5 Hz in the frequency range 0-140 Hz.

How can the Power Spectral Density of the combined digital signal that shown above can be estimated?


texts: the dominant frequency is overwritten in every iteration using matlab code

texts:
am trying to understand how to find a dominant frequency in a signal
i've a csvfile for an analog data wich detects the breathing rate for the human
i need to find the dominant frequency between 0.1hz and 2hz in this data for every 100 sample , so i wrote a matlab code to do so , after i run the code i find it gives me the same value in every iteration , so i think it must overwrite the value in every index,, here is the code

% clear the current variables in workspace
clear all;
close all;
clc;
filename = 'breathing.txt';
try
    T = csvread(filename);
catch
    warning('no such file');
end

% get the requested signals to process
signalx = T(:,15); % its length is 33167*1

%plot(signalx);

 % ----- apply median filter on data ------

windowSize = 5;
signalx = medfilt1(signalx, windowSize);
plot(signalx);

% -------- process the three signals -------

% initialize variables
samplingRate = 100;
movingwindowSize = 10*samplingRate;
incrementalWindow = 1*samplingRate;
% % % % tolerancefreq = 0.05;
dataLength = size(T,1);
dominantVector = zeros(floor(dataLength/movingwindowSize), 1);
count = 1;
% % % %
for i =1 : incrementalWindow : dataLength-movingwindowSize

    % ------ data segmentation -----

    currentwindow_sigx=signalx(i:i+movingwindowSize-1);

     % ----- apply hamming window -----

     hammingDatax=currentwindow_sigx.*hamming(length(currentwindow_sigx));
     plot(hammingDatax);

     %---- remove dc offset from data ----

     data_x = detrend(hammingDatax);

      % ------- apply fft ------------

     fftData_x = fft(data_x);
     plot(abs(fftData_x));

     % ------- find dominant frequency -----

     len = (length(fftData_x)/2)+1;
     datax = fftData_x(1:len);
     plot(abs(datax));

     fftLength = length(fftData_x);
     freq = (0:(samplingRate/fftLength) :(samplingRate/2));
     plot(freq,abs(datax)), grid('on'), title('dominant freq')

        freq =  freq(find(freq>=0.1 & freq<=2));

    [v,k] = max(datax);

      dominantFrequency = freq(k);
      dominantVector(count) = dominantFrequency;
     count = count +1;

end

so when i run the code it gives me a 322*1 vector with all values == .2000
am sure there is something wrong with the foor loop and the index but i can't find where te problem is.
and here is a photo of the workspace and sample of the overwritten values.


texts: Contradictory results from power spectral densities and autocorrelation in shot noise

texts:
I have a white noise generator circuit which involves a Zener diode meant to undergo breakdown across a resistor $R1$, and am trying to compare power spectral densities and autocorrelation of the output I get. I recorded the data on an oscilloscope with microsecond, milisecond, and half second timescales.

First, I tried looking at what sample rates that I chose from was going to best characterise the noise. I figured a lower sampling time would be better, sort of arbitrarily, but the noise histograms I got for lower sampling rates tended to look more Poissonian in nature (well, most of all Gaussian looking):

I guess I should take it that longer sampling rates tend to characterize the shot noise better?

As I say, Zener diode breakdown noise is shot noise, so I was ideally trying to see which combination of Zener diode voltage and $R1$ I could use to generate the "whitest" looking noise. That is, that could generate a flat power spectral density and most random-looking autocorrelation (all lags except for $1$ near $0$).

However, I have some contradicting graphs when I plot them. For instance, for a Zener diode voltage of $12 V$ and $R1 = 100 k \Omega$ my PSD looks like this (using plt.psd(data,Fs=1e6))

which is somewhat flat? However, the autocorrelation certainly disagrees with this being perfectly random (the axes are the canonical ones - lags and ACF):

which, if I am to assume the PSD shows flatness, contradicts it. Meanwhile, for a milisecond time between samples with the same Zener voltage and $R1$ resistance..

which again seems to contradict itself, but now in the other direction -- the near 0 autocorrelation supports the idea that noise is random, yet the PSD is not flat.

Meanwhile, with a half-second time between samples..

my data seems to be the most self-consistent here, if I can say that the PSD is more-or-less flat. Is this essentially telling me that the proper time between samples should've been one half-second? If not, what's going on here?

I'm very shaky to the idea that the faster sampling rate failed to capture the noise distribution as well as a slower one.


texts: How do I test stability of a MIMO system?

texts:
Let's say I have a system similar to two interconnected IIR filters described like this:

\begin{align}
x_1(t)&=a_{11} x_1(t-1)+a_{12} x_1(t-2) +a_{13} x_2(t-1) + a_{14} x_2(t-2)+y_1(t)\\
x_2(t)&=a_{21} x_2(t-1)+a_{22} x_2(t-2) +a_{23} x_1(t-1) + a_{24} x_1(t-2)+y_2(t)
\end{align}

where $x$ are the outputs and $y$ are the inputs.

How would I go about testing the stability of such a system?

I myself have pretty much no background in MIMO systems, so I'm not really sure where even to start.

For a SISO system, I'd start with poles and zeros and go on from there, but I'm not sure how such a graph might looks like in the MIMO case.


texts:
I would like some advice on which beamforming method to use from someone with experience in the field.

I have radar data from the ocean which has been transmitted by an omnidirectional transmitter and received by a linear array of 8 antennas with half wavelength spacing.

The signals come from all directions and I need to separate them using beamforming.

I started by using a simple phase shift beamformer which I wrote in Python but it wasn't filtering the signals from other directions well.

I've found the Capon beamforming technique to be an improvement. There's a toolbox on MATLAB which provides good functions for doing this.

However, my questions are:

What beamforming method would be best for this job?
What are the pros and cons of the different methods?

If any more details are needed, please ask.


texts: audio sound localization using machine learning

texts: I'm trying to localize a sound by using an array of microphones, the localization is done based on the sum and delay algorithm. I'm trying to to do the same in machine learning. I don't have much knowledge about machine learning but I think it can be done using linear regression using stochastic gradient ascend, I'm using a sine wave as input, I have to optimize the delay to and find out at which delay I get the max gain at the output. can anyone provide an example for linear regression using stochastic gradient ascend


texts:
In the paper "An Approach Based on Hybrid Genetic Algorithm Applied to
Image Denoising Problem", I encounter this fitness function that by minimizing it, a denoised image of a noisy image is obtained.

It is called Discretized Beltrami energy. According to the paper, $I$ is the image being evaluated, $I_0$ is noisy image, $\nabla I$ is total variation (TV) regularizing term, $\beta$ and $\lambda$ are balancing parameters. $\Omega$ is the set of all points of the image.

What is $\nabla I$? Determinant of $\nabla I$ is a constant, why is there a sum over $\Omega$?


texts: What is the frequency response function for an input through a differentiator and mixer?

texts:
Through the differentiator, the frequency response will be $j \omega X(j\omega)$, but what about through a mixer with $\sin(\omega_ct)?$

Will it be
$$\frac{1}{2j} j\omega\left[X(j(\omega-\omega_c))-X(j(\omega-\omega_c))\right]$$
$$\text{or}$$
$$\frac{1}{2j} j(\omega-\omega_c)\left[X(j(\omega-\omega_c))-X(j(\omega-\omega_c))\right]$$

I'm not sure if for a mixer you replace all $\omega = \omega-\omega_c$ or not.


texts: How to Smooth Gradient Estimates for Steepest Descent Optimization

texts:
In steepest descent
methods of minimizing a function $f(x), x \in \mathbb{R}^d$,

it's common to approximate the gradient by finite differences:
$\qquad\qquad \nabla f(x) \approx gradest( x; h ) \equiv
{ {f( x + h I ) \ - \ f(x)} \over h } $
then minimize along those directions:
$\qquad\qquad x_{k+1} \  \text{from min} \ f( x_k + \lambda \ gradest( x_k ) ) $

Can a signal-processing point of view suggest ways of smoothing $gradest()$ in this context ?
For example, one could view the "zig-zags" in the picture above as high frequencies in $gradest()$,
to be smoothed out over the last few ${x_k}$.
But that may be naive:
can one model, then smooth, zig-zags through non-uniformaly spaced points in 5d or 10d ?


texts: Calculation of autocorrelation - does calculation "loop"?

texts:
With a delay $l$, autocorrelation is defined as:
$$r_{xx}(l) = \sum_{n=-\infty}^{\infty}x(n)x(n-l) = \sum_{n=-\infty}^{\infty}x(n)x(n+l).$$

I want to calculate the autocorrelation of a signal
$$x(n) = [1,2,3,4,5]$$
with delay $l=1$. The sum starts of as
$$1*2+2*3+3*4+4*5..., $$
but what happens when I reach $n=5$? Do I start of in the beginning again, i.e is the autocorrelation
$$1*2+2*3+3*4+4*5 + 5*1$$
? Am I done there?


texts: How to detect the maximum resolvable spatial frequency of camera?

texts:
I am trying to calculate the minimum line pixel width that can be distinguished from noise as shown in the camera test chart in Figure 1 where the thinner lines on the left are getting more and more blurry and mixed with noise thus become unresolvable. I took the picture of the camera quality test chart using an LG V30 smartphone. The ideal image should look like Figure 2.

Does anyone have an idea about how to mathematically calculate the frequency or pixel size at which the lines are unresolvable or where lines becomes visible as noise?

I figured that the Fourier transform can be used (which is the bottom plot of Figure 1, but still struggling to know how).

Figure 1

Figure 2


texts:
I have the following aerial image (apparently taken from an airplane):

I also have a map of the area shown in the image:

On the image and the map I can mark the corresponding streets, street intersections and landmarks (I have shown only three here, but I can provide more points for precision).

I would like to skew the map to match the image (camera position, angle, etc).  How can I do this?

(I can program OpenCV and AForge.Net, and can use Gimp and Blender, all on Windows.)


texts:
It is documented that 'one' of the units of the Fourier Transform [of $x(t)$ volt] is volt per Hz. That is $X(\omega)$ components will have units of volt per Hz, where $\omega$ is the angular frequency (along the $x$-axis).

For example, the Fourier Transform of $\cos(\omega_0 t)$ volt has two components - one at $\omega = \omega_0$ and one at $\omega = -\omega_0$, each assigned a 'value' of $\pi$. So it is assumed that the values are '$\pi$ volt per Hz'.

Now, when the plot is adjusted to become a function of cyclic frequency $f$ (in Hz), which means plotting $X(f)$ versus $f$, then the Fourier Transform components of $\cos( 2\pi f_0 t)$ become '1/2 or 0.5' instead of '$\pi$'.

My question is --- why isn't the '$\pi$ volt per hz' preserved when translating between $X(\omega)$ and $X(f)$?

According to tables showing Fourier Transform pairs, the Fourier Transform of  $\cos( 2\pi f_0 t)$ volt has two components -- each with a value of '1/2 or 0.5'. For these values, does anyone know what are their units? I had assumed that the values for each component would have been 'preserved' (to be $\pi$ volt per Hz), regardless of whether we do the plot as a function of '$\omega_0$' rad/s or a function of '$f_0$' Hz. Thanks all.


texts:
Relatively new to DSP, I've came across the approach for phase-estimation of PSK-Signals x(k):

My question is:

What do we get by exponentiating the Baseband-Signal?


texts: Understanding Integration Scale in Multi-Scale Harris Corner Detector

texts:
I would like to improve a simple Harris corner detector by adding more scales.

I don't need scale-invariant interest points, just to detect more points by processing more scales.

The Wikipedia says there are two sigmas involved: a derivative scale $\sigma_{D}$ and an integration scale $\sigma_{I}$.

So after computing derivative images on a given level, I smooth them with Gaussian filter of size $\sigma_{D}$, then for each location $(x,y)$ in the image I have the matrix:

$H(x,y)=\begin{pmatrix}I_{x}^{2} & I_{x}I_{y} \\ I_{x}I_{y} & I_{y}^{2} \end{pmatrix}$

where $I_{x}$ and $I_{y}$ are the smoothed images (with $\sigma_{D})$ on the respective location.

How do I apply the integration scale $\sigma_{I}$ now?

Do I have to store separate images $I_{x}^{2}$, $I_{y}^{2}$ and $I_{x}I_{y}$ and smooth these with $\sigma_{I}$ before computing the matrix and Harris corner response from it?


texts:
Can someone help me:

I have a signal which consist of three frequencies; for example, $f_1=800\,\text{Hz}$, $f_2=1050\,\text{Hz}$, $f_3=1600\,\text{Hz}$. How can i compute the fundamental frequency?

Also, is there a method called "frequency histogram method" which is used to compute fundamental frequency? Can someone explain it to me?

Thanks.


texts: Significance of $ \lambda $ in Basis Pursuit

texts:
In basis Pursuit, L1 minimization is done to perform compressed sensing. In the literature there is a $ \lambda $ parameter used as a regularizer.

What is its significance?


texts: what signed 16 bits little endian mean in pcm audio buffer?

texts: I'm trying to read data, I want this data in numbers (int double float something i can operate with in order to make some dsp) from a buffer obtained from a pcm configured as S16_LE, (signed 16 bits little endian) but i don't understand how information is organized inside the buffer. I think 16 bits means that signal values are between -32767 and 32768 but nothing more..


texts:
This question stems from a point of confusion that I still have about the causality, linearity, and time-invariance in LCCDEs. I wanted to use the capacitor as an example.

Consider a capacitor with capacitance $C$. Taking the current $i(t)$ to be the input to the system and the voltage $v(t)$ to be the output we have $$i(t) = C \frac{\mathrm{d} v(t)}{\mathrm{d}t}$$

This differential equation can be solved to obtain
 $$v(t) = v(t_0) + \frac{1}{C} \int_{t_0}^{t} i(\tau) \mathrm{d} \tau$$

My first question is: isn't this mathematically valid for all $t$? In other words, does this give us the response for all $t$ or is it only valid for $t > t_0$? If it is valid for all $t$, including the $t < t_0$ case, doesn't this make the system non-causal since it anticipates future input and output values? Are we allowed to integrate backwards in time?

My second question relates to the assertion that a for the LCCDE to describe a linear system, the initial conditions must be zero. Suppose $t_0 = 0$ such that  $$v(t) = v(0) + \frac{1}{C} \int_{0}^{t} i(\tau) \mathrm{d} \tau$$

With $v(0) = 0$ the system is linear. But the choice of $t_0 = 0$ is arbitrary, since for example
$$v(t) = v(0) + \frac{1}{C} \int_{0}^{t} i(\tau) \mathrm{d} \tau = v(2) + \frac{1}{C} \int_{2}^{t} i(\tau) \mathrm{d} \tau$$

Why shouldn't we require that $v(2) = 0$ as well for that matter? What am I missing here? Thank you in advance.


texts: How to make semilogx to get a plot gain against log frequency with several values? [beginner]

texts:
Gain: 20 log (Vout/Vin) with Vin = 1.0. I've been searching for gain values: a. -156.48 dB b. -80 dB c. -53.98 dB d. -13.98 dB e. 3.52 dB f. 10.32 dB g. 13 dB h. 13.93 dB i. 13.96 dB j. 13.98 dB

I am still at the very start on this semester to learn about matlab. How to make semilogx to get a plot gain against log frequency with several values? I've been searching matlab tutorial but I still got no clues. If anyone could give some explainations or examples relate with this question, please tell me. I want to learn about it. Thank you.


texts: why wavelet coefficients are zeros after decomposition

texts:
I am learning wavelet theory for image processing. To understand the theory, I write one Matlab program to decompose one black-white image. The program is as follows

 Image = zeros(256, 256, 'uint8');
 Image(101:200, 101:200) = 255;
 figure; imshow(Image);
 [cA1, cH1, cV1, cD1] = dwt2(Image, 'db1');
 Image1 = [cA1, cH1; cV1, cD1];
 figure; imshow(Image1, []);
 [cA1, cH1, cV1, cD1] = dwt2(Image, 'db2');
 Image1 = [cA1, cH1; cV1, cD1];
 figure; imshow(Image1, []);

The first decomposition using the argument db1 produces zeros for all wavelet coefficients. The black-white image has the transition from 0 to 255 along horizontal and vertical directions and should have high-frequency component. Why are zero wavelet coefficients generated?
If I change the argument from db1 to db2, the result will show horizontal and vertical lines in the subbands.

Edit:
I decompose the image further and at level 3 edge-like structure can be seen on the subband. Wavelet is one kind of filter bank I read from one book. It means the frequency in the image can only be seen in the higher frequency band. From the below link, it told how to get frequency response of Daubechies filters.

http://www.mathworks.com/matlabcentral/answers/44879-daubechies-filters-frequency-response

I don't understand how is the frequency for images defined. From one book, it says that "Spatial frequencies are measures of the number of oscillations (cycles) from dark to light per unit of length of the image". How can I know the frequencies in the input image.


texts:
In standard of 3GPP LTE physical layer(3GPP TS 36.211), the FFT size for baseband signal generation for both Dl and UL is 2048 and CP size is also determined. However there is a table in several resources said we have different FFT and CP size. E.g.

http://lteuniversity.com/ask_the_expert/f/59/t/405.aspx

How is it possible?


texts:
I have a sound file that needs filtering which looks like this:

Hence I think that I have noise spaced equally at 832Hz. i tried to use the iircomb filter taken from mathworks website and i managed to create this:

clear
fs = 44100;
fo = fs/53;
q = 35;
bw = (fo/(fs/2))/q;
[b,a] = iircomb(fs/fo,bw,'notch'); % Note type flag 'notch'
fvtool(b,a);
[f,fs] = audioread('sample3b_va18535.wav');
fnew=filter(b,a,f)

N = size(f,1);
df = fs / N;
w = (-(N/2):(N/2)-1)*df;
y = fft(fnew(:,1), N) / N; % For normalizing, but not needed for our analysis
y2 = fftshift(y);
figure(2);
plot(w,abs(y2));
p = audioplayer(fnew,fs);
p.play;

Now, the filter works pretty good but there is still a little bit of noise after filtering which can be heard. Should I just live with it, or perhaps my filter could be improved? Could someone with more experience work on it for a minute or two? The sound file is only 3 seconds long.

Here is the audio file if anyone needs it: SOUND FILE

Thank you for help.


texts: Inverse Discrete-Time Fourier Transform of $X(Ω)=jΩ$

texts: I am trying to solve it by using the properties but I can’t seem to find the same solution as on my book.


texts: Interleaving over OFDM symbols in same transmission

texts:
In this question I have asked what are the benefits from interleaving between two codewords over ONE OFDM symbol in packet based transmission made up of more than one OFDM symbol

Purpose of interleaving in OFDM systems

The answer was satisfying. I have another question now and I am looking for your opinion on this.

If I decide to do interleaving outside of a single OFDM symbol, (say for example in packet based transmission) then is there any benefit?

Thanks


texts: The complexity of such function run in Matlab

texts:
The below function is representing an algorithm, so how can I get its complexity?  I don't mean the time of running by using the tic .. toc, I mean how many operation (Additions and multiplications) are performed in this loop.

for times=1:m;
    for col=1:N;
        product(col)=abs(T(:,col)'*r_n);
    end
    [val,pos]=max(product);
    Aug_t=[Aug_t,T(:,pos)];
    T(:,pos)=zeros(M,1);
    aug_y=(Aug_t'*Aug_t)^(-1)*Aug_t'*Yy;
    r_n=Yy-Aug_t*aug_y;
    pos_array(times)=pos;
end

Size of parameters, m = 256 , N = 256, T= [256,256] and M = [256,1]


texts:
If we know that a filter (or a system) has two poles at 20 GHz and one zero at 15 GHz, then how do you write the transfer function $H(s)$ for such a system?

I am wondering why sometimes the poles and zeros are given in Hz and other times as complex numbers?


texts: Why call it discrete time signals as difference equations?

texts:
There are lots of representations of discrete time signals. You represent it in block diagram or as difference equations.

We call it block diagrams as there are elementary structures in the diagrams.

Why do we call the equation part as difference equation of a signal?
$y[n] = x[n] + x[n-1]$

There is nothing getting subtracted in the above example representation, Then why call it difference equation??


texts: Identify the system from input and output

texts: I am using Excel. If I have an input and an output, and I know that a recursive band-pass filter produced the output, how can I can find the coefficients of that filter?


texts:
I am researching downscaling methods for image downscaling. I've both read and found area-averaging to work well in Python using CV2, but the settings in MATLAB allow for it's bicubic interpolation downscaling to work better.
Are there any resources comparing image interpolation techniques? especially explaining the theory? (I can't seem to find any books or journal articles that cover what I need).
Thanks


texts: Differences between Python and MATLAB filtfilt function

texts:
I'm trying to port some MATLAB code to Python and am running into some strange behaviour. I am implementing a 5th order Butterworth bandpass filter. The sampling rate is 30 Hz.

Running MATLAB R2012b, Spyder 2.2.0 with Python 2.7, SciPy 0.12.0 on Windows 7 x64.

In MATLAB:

[b,a] = butter(5, [0.75*2/30, 5.0*2/30], 'bandpass');
y = filtfilt(b, a, input_signal)

This is the raw signal:

and the filtered signal:

and its power spectrum:

which makes sense since the normalized bandpass frequencies are 0.05 - 0.33.

I found before that SciPy's butter function does not give the same coefficients as MATLAB so I exported the filter coefficients from MATLAB to Python using hdf5write (see here: https://stackoverflow.com/questions/7117797/export-matlab-variable-to-text-for-python-usage)

In Python:

y = signal.filtfilt(b, a, input_signal, padtype = None)

and the output is:

and its power spectrum:

I used padtype = None because by default it is padtype = 'odd'. However, I've tried all the different padding options and they all look more or less the same.

I'm not entirely sure what's going wrong...any help would be greatly appreciated.

EDIT: Added graph for padtype = "odd" and b and a filter coefficients used

Signal with padtype = "odd":

Seems like there is some underlying signal so I don't think it's the impulse response although that large transient at the beginning is strange.

Python filter coefficients (5th order Butterworth filter):

passband = [0.75*2/30, 5.0*2/30]
b, a = scipy.signal.butter(5, passband, 'bandpass')

b, a are arrays of type float64.

b = array([  5.49209388e-03,   0.00000000e+00,  -2.74604694e-02,
    -1.97776791e-17,   5.49209388e-02,   2.47220989e-17,
    -5.49209388e-02,  -1.97776791e-17,   2.74604694e-02,
     0.00000000e+00,  -5.49209388e-03])
a = array([  1.        ,  -6.52098852,  19.50384534, -35.47189804,
    43.65758795, -38.07760914,  23.83047021, -10.55670367,
     3.16710078,  -0.58122912,   0.04945954])

MATLAB filter coefficients

b =   0.0055, 0, -0.0275, 0, 0.0549, 0, -0.0549, 0, 0.0275, 0, -0.0055
a =   1.0000, -6.5210, 19.5038, -35.4719, 43.6576, -38.0776, 23.8305,
      -10.5567, 3.1671, -0.5812, 0.0495

Turns out the coefficients are about the same, I believe the issues I had before were due to the frequency being much higher (coefficients must be generated differently in MATLAB and Python).

The input signal is a list of float32, although I've tried converting to array with numpy.array and the result is the same.

EDIT: More information about padtypes and Python vs. MATLAB

Python SciPy's filtfilt function includes a parameter called padtype which indicates the type of padding extended on both sides of the signal. This padding serves to reduce transients. Odd and even are descriptors of the type of symmetry these extensions have with the endpoints (http://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.filtfilt.html).

Some illustrations of the difference in padtypes:

padtype = "even"

padtype = "odd"

padtype = None

Based on the results, it appears that MATLAB uses odd padding (which is also the default for Python):


texts: Are the Daubechies 4 wavelet function and scaling function dimensionless?

texts: If the dimension of the independent variable $x$ of $f(x)$ is length, i.e. $[x]=L$, then what is the dimension of Daubechies 4 (hereafter D4) wavelet and scaling function?


texts: How to count the number of threads in a bolt or nut using image processing?

texts:
How can I count the number of threads in a Nut with image processing?

I know that I need to find the place of threads and count them. Ok, But how can I find the place of threads?

This is a sample image:

Update
As I know I need to do these steps:

Making image ready. For example az @A_A sayed, with Thresholding
Remove noise and just keep threads in the image
Count threads. maybe using signal processing

Level 1 is solved. What should I do for 2 and 3?

For example this is what I did for a bolt. I got a sample of a vertical line in the middle of the nut then plotted the data. But I don't know how to remove additional parts and count up and downs.

Thanks


texts: Wiener Filters - Obtaining / Estimating Image and Noise Spectrum

texts:
We all know the famous Wiener Filter,

$$G(\omega_1 , \omega_2)=\dfrac{H^*(\omega_1 , \omega_2)}{|H(\omega_1 , \omega_2)|^2 + \dfrac{N(\omega_1 , \omega_2)}{S(\omega_1 , \omega_2)}}$$

How to we use it practically? Meaning, given $y = s * h + n$ where $s$ is the original image, $h$ is a LSI operator, and $n$ is some noise (independent of $s$), how to calulate (or estimate) $N(\omega_1 , \omega_2)$ and $S(\omega_1 , \omega_2)$ in order to reconstruct $s$ from $y$?


texts: TypeError: primitive_connect() when executing in gnu radio

texts:
I'm building the frequency analyzer for RTL-SDR. But having problem in building and executing code.
QStandardPaths: XDG_RUNTIME_DIR not set, defaulting to '/tmp/runtime-root'
gr-osmosdr 0.2.0.0 (0.2.0) gnuradio 3.8.1.0
built-in source types: file osmosdr fcd rtl rtl_tcp uhd miri hackrf bladerf rfspace airspy airspyhf soapy redpitaya freesrp
[INFO] [UHD] linux; GNU C++ version 9.2.1 20200304; Boost_107100; UHD_3.15.0.0-2build5
Found Rafael Micro R820T tuner
Using device #0 Realtek RTL2838UHIDIR SN: 00000001
Found Rafael Micro R820T tuner
[R82XX] PLL not locked!
[R82XX] PLL not locked!
Traceback (most recent call last):
  File "a.py", line 209, in <module>
    main()
  File "a.py", line 187, in main
    tb = top_block_cls()
  File "a.py", line 144, in __init__
    self.connect((self.rtlsdr_source_0, 0), (self.qtgui_freq_sink_x_0, 0))
  File "/usr/lib/python3/dist-packages/gnuradio/gr/hier_block2.py", line 37, in wrapped
    func(self, src, src_port, dst, dst_port)
  File "/usr/lib/python3/dist-packages/gnuradio/gr/hier_block2.py", line 100, in connect
    self.primitive_connect(*args)
TypeError: primitive_connect(): incompatible function arguments. The following argument types are supported:
    1. (self: gnuradio.gr.gr_python.hier_block2_pb, block: gnuradio.gr.gr_python.basic_block) -> None
    2. (self: gnuradio.gr.gr_python.hier_block2_pb, src: gnuradio.gr.gr_python.basic_block, src_port: int, dst: gnuradio.gr.gr_python.basic_block, dst_port: int) -> None

Invoked with: <gnuradio.gr.gr_python.top_block_pb object at 0x7fe58fbaf630>, <Swig Object of type 'gr::basic_block_sptr *' at 0x7fe58bd46540>, 0, <gnuradio.gr.gr_python.basic_block object at 0x7fe58e3915b0>, 0
swig/python detected a memory leak of type 'gr::basic_block_sptr *', no destructor found.
```


texts:
Signal processing is Elsevier, so it's not open access, and IEEE Transactions on Signal Processing is also not open access. I was wondering if there's any decent journal on par with those that's open access.

I'm trying to submit a paper on discrete time nonlinear filtering.

(there are no tags relating to publication/journal/conference please feel free to add tags or suggest different SE)


texts: Using an IMU for gravimetry and calculation of gravity anomalies

texts:
I am currently working on a project which seeks to utilize a tri-axial IMU to measure gravity anomalies caused by density variations in planetary crusts.

The instrument will be attached to a Helikite balloon at a maximum altitude of approximately 30 meters and towed by a rover. I understand the errors associated with IMU data, however, I am having some difficulty understanding how to implement the different corrections to the data.

Is a rotation matrix/quaternions required in order to shift the gravity vectors into the earth/navigational frame?

Also, any advice on filter techniques to help account for drift (e.g. Kalman filter) would be immensely helpful. My field of knowledge is earth/planetary science so I'm not too familiar with these kinds of techniques.


texts: Channel coder: Does encoder help? Input signal distortion:

texts:
I am studying the channel coding theory. The channel coder adds redundant information to the input signal, that helps a decoder to reconstruct the signal.
BUT: Can encoder help to eliminate effects such as Doppler, offset...?
Edit 1
when we are talking about improving encoder in practice, what does it mean? What should be improved?
Edit 2:
Why has an interleaver complicated implementation in practice? How to reduce complexity?
Edit 3
S/P or P/S converters and interleaver, do they make the same task?


texts: Loop bandwidth for symbol timing recovery

texts:
I am trying to design the Symbol timing recovery (STR) block of DVBS Receiver.

How should I chose the loop bandwidth and damping ratio for determining the loop filter coefficients $k_1$ and $k_2$?

In some literature they are given as $B_nT_s= 0.005$, $\zeta = \frac{1}{\sqrt 2}$, while in others it is given as $B_nT_s= 0.05$, $\zeta = 1$.

How should I proceed for choosing the $B_n$ and $\zeta$?
Should $B_nT_s$ factor change as we change $T_s$ (say from $0.1$ to $10^{-5}$ sec)?


texts: Extended Kalman Filter in mechanics, electronics and hydraulics?

texts:
Extended Kalman Filter is most used in GPS and navigation systems. But how much do I gain to switch to Extended Kalman Filter to the linjear Kalman Filter if I do LQGI controller for the industry?

All physics are nonlinear like hydraulics and mechanics.
What do you prefer?


texts: What type of low-pass filter to use for gyroscope filtering in control application?

texts:
I am using a gyroscope to generate a feedback signal that drives a motor to resist rotation, just as in a camera gimbal.

The gyroscope is sampling at 32 KHz with a bandwidth of 8 KHz, and my control loop is running at 1 KHz. Hence, I need to remove frequencies above 500 Hz. At the same time, I need to introduce as little delay as possible in order to be able to use a high control gain. I've attached the FFT of the sensor noise (source).

I am a little overwhelmed with the amount of choice available when it comes to filtering this signal. What filter properties should I be looking for? Are there any particularly suitable filters? This is implemented on an FPGA, so computational cost is not a constraint.


texts:
I'm trying to get a handle on the importance of the error feedback term in noise shaping operation in typical audio dither algorithms.

I'm thinking in terms of four signals.  The original signal, the output (quantised) signal, the dither signal, and the error signal (the difference between output and original).

A simple approach would be to generate white noise and apply a filter to it to generate a coloured dither signal, add this to the original signal, quantise to the output signal, and disregard the error signal.  It appears to be a fairly common approach to use white (unfiltered) or violet (differentiated) noise functions in this configuration.

Another approach I see implemented is to take the error signal and filter that and add white noise to define the dither signal, then add it to the original, quantise, and calculate the next error term for the filter.  This implies a tight loop for the data flow -- we cannot compute a new error term without resolving the previous one -- making it hard to use SIMD efficiently.

From those I can infer a third option -- forgo the noise altogether and simply filter and reintegrate the error signal directly.  I suppose this comes with some risk of modulation artefacts related to but distinct from the artefacts introduced by quantisation without any dither or shaping at all; however, I suspect this might be how it's done in some high-frequency 1-bit cases.

The point I'm unclear on is that complex filters which attempt to model the threshold of human hearing appear to come under the title of 'noise shaping', and this implies (supported by implementations I have seen) use of the error feedback term.

Is it a flawed or misguided approach to forgo the error feedback but still use the complex filter on white noise for the dither function?  Is the error signal genuinely offering something important here?


texts: From transfer function to differential equation

texts:
I have the below detailed solution (boxed in blue) that I don't understand completely:

I can reconstitute the differential equation from:
$$ (1+Ts) X(s) = K_v U(s) $$
$$ x(t) + T\dot x(t) = K_v u(t) $$
but I don't get how do they get $x(t)$ from the initial conditions?
I realize that my question is more about differential equations but it is something I obviously missed before...


texts:
I know I can write the frequency response of a Generalized linear-phase system as:
$$H(e^{j\omega}) = A(\omega) e^{-j\left(\alpha \omega - \beta\right)}$$
where $A(\omega)$ is real.

I need to prove that $A(\omega)$ is real and I can't find anywhere the demonstration.


texts: i just wondering, why we always transmit a real signal but when we deal with a baseband signal we use a complex signal ? are this is related to up and down conversion of the signal (Because complex signal doesn't have symmetric properties when using FFT)


texts: What is the role of a LPF in oversampling?

texts: Lets assume that a signal $x(n)$ is up-sampled by adding 1 zero between two adjacent samples to form a signal $y(n)$. How does a digital LPF give an oversampled version $z(n)$, with more data points? (As show in the image below) I know that the question might seem stupid but the books that I am using don't really explain that part.


texts: Why Carrier signal can't have frequency less than Message signal

texts:
I know that why we use Carrier of high frequency to send message signal over a long distance, when the frequency is low , energy will be obviously low. To increase the energy of the signal we need to increase the frequency. This is achieved by multiplying the message signal with the carrier signal (with high frequency).

suppose,

Carrier signal frequency = 2800KHz
message signal frequency = 3KHz

the two generated sidebands will be ,

2800 + 3 = 2803 KHz
2800 - 3 = 2797 KHz

and eventually, the bandwidth of the signal is,

BW = 2803 - 2797 = 6KHz

this is when the carrier frequency is higher , and we noticed that bandwidth is just twice of the highest frequency of modulating signal.
and the signal will be easily demodulated at the receiver end.

but what will happen if it is reversed ?

i.e

Message signal frequency = 2800KHz
Carrie signal frequency = 3KHz

please explain would happen here ?


texts:
We know that a discrete-time system with a (Z-transform) transfer function that has a pole of magnitude 1 (i.e. $|z|=1$ is a pole of the transfer function) is marginally stable if the pole at $z=1$ is a single pole and unstable on every RoC otherwise.

My question is: Suppose a transfer function that has two simple poles at $z=1, z=-1$. Is it marginally stable or not? And, therefore, what is the most general form of the stability theorem based on poles of the (Z-)transfer function?


texts: Why can't we just make all wireless networks use integer multiples of base frequency?

texts: I always wondered why transmission capacity depends on bandwidth. For example, let us say that there is an isolated island. In this island, people decide that all wireless networks use frequencies that are positive integer multiples of one base frequency. Then when we do Fourier transform, we can just use least common multiples of all used frequencies and data in each frequency will exactly be determined except background noises. Or am I wrong here? In this sense, bandwidth seems just to be useless....


texts: Find the frequency response of the system

texts: How can I find the frequency response of the following linear system


texts: What is the exact meaning of the output of the Discrete Fourier Transform

texts: I'm fairly new to the subject, but so far my understanding that this would be a transform you could use to go from a discrete set of data, say [1, 0, 1, 2] to a continuous sinusoidal function in the frequency domain, if I have fundamentally got this wrong maybe this is where my confusion comes from but once you have applied the DFT to this signal, and then find its magnitude. What is the meaning of the outputs, are they coefficients of a sinusoidal function? Or is the just the representation of that data as an amplitude in accordance to its frequency?


texts: Fourier Transforms and Series for the NON mathematically inclined.

texts: This would most likely be the opposite of this question ( Mathematically inclined Signal and Systems/Signal Processing book? ) I figured I'd ask here if there are any good books that while, necessarily, focus on the math are more oriented to what these concepts can DO instead of the mathematics behind it, so I guess more how it changes the graphing or the function etc instead of the functions behind them.


texts: What's the number of data points in a Welch-based PSD?

texts:
In Welch periodogram Power Spectral Density estimate, we devide the N-long signal into K segments, each of length L with overlapping D, such that N=L+(K-1)*D.

In this paper, just above eqn (7), the authors mention that the step of performing the DFT is to be done to each (windowed) segment (of length L) such that the length of the resulting DFT is L/2. Finally all DFT's will be averaged and normalized.

At the end, you get a PSD of a max. length L/2. Why is that? I would expect the length of the resulting PSD to be related to N, not to L...


texts: Affine Transformation for matching 2D point sets

texts:
I have two sets of 18 corresponding 2D points. These coordinates represent the main joint points of a human body, describing the "human pose".

I want to decide if these two poses are similar/have the same shape. In order to do so, I find the affine transformation (matrix) that transforms the second pose to the first pose (using least-squares). This gives me the image of the second pose pictured on the first one.

Now, my question is, what is the best way to conclude if the image of the second pose is similar to the first one? I'm thinking about using a distance metric like mean square or SAD, but I'm not sure. + What about normalizing/standardization?


texts: How do I make a Gaussian pulse that works on different wavelengths in MATLAB?

texts:
I want to make this MATLAB code to plot a Gaussian pulse in THz, then use FFT on it. This is the MATLAB code I used to get the FFT of time domain Gaussian signal at 800nm wavelength:

lambda= 800e-9;%Wavelength
c=3*10.^8;%lLight speed
fc = c/lambda; % Actual Frequency of light [THz]
fsamp=fc*10; %sampling frequency
fs=1/fsamp; %Unit Time in [fs]
Ls=200; %Length of signal
sigma=8e-15;% Pulse duration
ts=(0:Ls-1)*fs; %time base
t0 = max(ts)/2; % Used to centering the pulse
gauss = (exp(-2*log(2)*(ts-t0).^2/(sigma)^2)).*cos(-2*pi*fc*(ts-t0));
subplot(2,1,1)
plot(ts,real(gauss),'b');
title('Gaussian Pulse in Time Domain');
xlabel('Time(s)');
ylabel('Amplitude');
grid on;
NFFT = 2^nextpow2(Ls);
X = fft(gauss,NFFT)/Ls;
freq = 0.5*fsamp*linspace(0,1,NFFT/2+1); % (full range) Frequency Vector
subplot(2,1,2)
plot(freq,abs(X(1:NFFT/2+1)))
title('FFT of the Gaussian Pulse');
xlabel('Frequency (Hz)')
ylabel('Magnitude |X(f)|');
grid on;

Which will give:

However, when I change the wavelength to 30um, it will give me a different shape:

I know this is due the pulse duration (sigma), or maybe entirely wrong code. How do I make it retain its shape whenever I change the wavelength, without continuously changing the pulse duration?


texts: Background subtraction for object tracking

texts: i am working on object tracking. i want to use background subtraction. i have tried to find out papers on background subtraction but i am really confused to use which one. I have already done tracking using just gray scale subtraction and thresholding. but now i want to model the background and then use it for tracking. i do not know much about background modelling.So, if anyone can suggest me some techniques and good papers on background subtraction and modelling, i will be greatly thankful.


texts: 3D point cloud to 2.5D graycale image (metric) [Halcon?]

texts:
I have a problem ragarding transferring a point cloud to a 2.5D metric grayscale image (gray value gives depth information). I start with a three channel image of 32bit floats with each pixel representing the three coordinates in mm (metric).

Halcon is the preferred library, but not necessary.

I hope that someone can help me.

Greetings from Stuttgart, Germany

Max


texts: Are there any direct relation formula between BER and Transmitter gain. The higher the transmitter gain is, the better the BER performance is. In here, I want to plot graph for the BER based on different transmitter gains and compare the MATLAB simulation vs experimental results.


texts: Making Unscented Kalman Filter Robust for Nonlinear Parameter Estimation Problems

texts:
So I have built code for an Unscented Kalman filter that can take any specified state and measurement dynamics. I have tested it on various linear problems and it works well, as expected. The main concerns I have is applying this filter to nonlinear parameter estimation problems.

When I have applied it to problems where I give a pretty good guess, it tends to do fine. However, more often when I don't have a good guess and make the covariance matrix elements larger, it tends to lose its positive definiteness and in turn fails to work.

Does anyone have tips on helping avoid these problems other than trial and error tuning of the various covariance matrices?

I have read papers about using the EKF/UKF for parameter estimation of weights in Neural Networks, for example. However, whenever I try to do this, I lose positive definiteness and it's frustrating because I struggle to reproduce the stellar results I see in the literature with respect to using the UKF for machine learning. Due to my failed attempts, I am questioning if it's really as useful as the literature appears to imply.

Note that I tend to formulate the problem like so for the parameter estimation problems I am looking at:

$$
\begin{align}
\textbf{w}_{i+1} &= \textbf{w}_{i} + \mathbf{\eta}_w\\
\textbf{e}_{i} &= \text{H}(\textbf{w}_{i}) + \mathbf{\eta}_e
\end{align}
$$

where $\textbf{w}$ is the parameters being estimated and $\textbf{e}$ is some error generated by the parameters, where the expected value for $\textbf{e}$ is always $\textbf{0}$ (since I'm trying to find $\textbf{w}$ such that the error is as close to 0 as possible). $\text{H}(\cdot)$ also tends to represent a cost function that the underlying model is scored against as a function of the current parameters $\textbf{w}$.


texts: Computing response of LTI system in terms of its step response

texts:
I have solved a problem, Kindly help me in determing if I solved it correctly.
I will post the question and my own working below.

Problem

My working


texts: Signal's Fundamental Frequency is different from Plotted Signal

texts:
I've been attempting to plot the following function using MATLAB:
$$
x(k)=\sum_{n=11}^{50} \sqrt{n} \sin (2n\pi k) +\sum_{n=1}^{40}\sqrt[3]{n} \sin (3n\pi k)
$$

Note that $k$ is a continuous variable.

The problem I am facing with MATLAB is when I plot it, I get a signal which has a fundamental period of 2 seconds, and not 1.5 seconds.

Have a look below:

I've tried a higher sampling rate:

So, what is going on wrong? Why isn't the fundamental period 1.5Hz?


texts: What would be the effect of repeating samples in time domain on frequency domain?

texts: Consider the operation of up-sampling a signal n times by repeating each sample n times in the time domain. What would the frequency domain response of the resultant signal be like given the frequency response of the original signal?


texts: How to setup filters to catch 4FSK deviations

texts:
I'm trying to demodulate a radio signal from Motorol DP4801 (MOTOTRBO) that follows DMR specification.

I'm receiving the signal by using HackRF with the rate sample 1 MHz. The radio uses frequency 406MHz. According to the standard DMR uses 4FSK modulation with the following deviations:

Symbol +3 - +1,944 kHz

Symbol +1 - +0,648 kHz

Symbol -1 - -0,648 kHz

Symbol -3 - -1,944 kHz

The RF carrier bandwidth is 12,5 kHz.

Because after the sampling the radio signal I've lost the original signal (406 MHz) I can't understand how I should set up bandpass filters to catch these deviations. I need it because I'm going to apply non-coherent FSK demodulation but as in the image but for 4 symbols.


texts:
I am trying to wrap my head around the original wavelet shrinkage papers of Donoho and others from the mid-90s. There are a couple of curious things that confuse me, and I'd appreciate some input:

The universal threshold is often (and also in secondary literature) given as $\sigma\sqrt{2\ln n}$, but sometimes I see $\frac{\sigma}{\sqrt{n}}\sqrt{2\ln n}$, e.g. in "Wavelet Shrinkage: Asymptopia?". The latter seems utterly wrong...
The minimax threshold $\lambda^\star_n$ as defined in "Ideal spatial adaptation by wavelet shrinkage" has the property that $(\lambda^\star_n)^2 = 2\ln(n+1) - 4\ln(\ln(n +1)) - \ln 2\pi + o(1)$. I've tried to use the square root of this to derive values for $\lambda^\star_n$  (ignoring the $o(1)$, guess that's what's causing my problem?), however, they are very different from the ones in table 2. Is there an explicit formula or algorithm for $\lambda^\star_n$?


texts: Role of Kalman filter in nonlinear dynamics

texts: I am quite interested to know the significance of using kalman filter,unscented kalman filter and extended kalman filter in chaos synchronization when infact the very basics of chaos synchronization(physics perspective) never mentions about using kalman filter!However, these papers KF1 and alike use kalman filter for synchronization. There are various methods of chaos synchronization viz generalized, phase, projective, adaptive Boccaletti_Kurths et al and wikipedia considered to be the bible in thi field, do not mention about using kalman filtering and alike methods. Would be obliged if the role of this technique is explained since I think I am missing on important theory. Thank you


texts:
I want to ask Question about the Fourier series in continuous time domain while reading a book signals and systems Alan Oppenheim. I have confusion in understanding the statement on page 189 of its 2nd edition.
To derive alternative form of Fourier series, we first rearrange the summation in Eq.3.25 as
$$x(t) = \sum^{+\infty}_{k\ =\ 1} [a_k e^{jk\omega_0t}+a_{-k} e^{-jk\omega_0t}]$$

where Eq. 3.25$$x(t) = \sum^{+\infty}_{k\ =\ -\infty} a_k e^{jk\omega_0t}$$
I want to know two things
1) How summation from -infinity to +infinity changes to 1 to +infinity

2) How we get this term $$a_{-k} e^{-jk\omega_0t}$$ in the equation.


texts: Why is $e^{i2\pi Nf_snT}=1$ in snippet from book?


texts:
I am studying control systems, and have encountered the definition of a slow zero. I am searching on internet and in books, but I don't understand the meaning this. I know that if a zero is too slow, it introduces a resonance peak and an overshoot in the step response, but what is a slow zero?
I thought it was a zero close to the imaginary axis in the root locus, so a zero at low frequencies, but I am not sure, since I cannot find a definition.
I know that a pole that is close from the imaginary axis settles quicker than a pole far from the imaginary axis. But I did some simulations in Matlab, and using a lead compensator, I have seen that if I decrease frequency of the zero in the lead compensator, the overshoot in the step response increases.

Here, the lead compensators I have used are as below.
lead = (1+10*0.05*s)/((1+0.05*s));          %red line
lead_2 = (1+30*0.083*s)/((1+0.083*s));      %blu line
lead_3 = (1+50*0.083*s)/((1+0.083*s));      %green line

Can somebody help me?


texts: Make a Movie from Multiple Images in MATLAB

texts:
I have 10 images in the folder say "D:\images\" named as frame1,frame2,...,frame10. I want to make a movie from the color images with frame rate as per choice.Also i want to save created video file in .mp4 format at the location "D:\video\".

Is it possible to do so in MATLAB?how?


texts: Computing Discrete Convolution in terms of unit step function

texts:
Consider :
$$
x[n]:=\begin{cases}
1&\text{if $3\leq n\leq 8$}\\
0&\text{if otherwise}
\end{cases}
\quad\text{and}\quad
h[n]:=\begin{cases}
1&\text{if $4\leq n \leq 15$}\\
0&\text{if otherwise}
\end{cases}
$$
I wish to compute the convolution $\displaystyle y[n]=(x*h)[n]=\sum_{k=-\infty}^{\infty}x[k]h[n-k]$. Because I dislike the method of graphing, I wanted to take another approach using unit step function :
$$
x[n]=u[n-3]-u[n-8]
\quad
\text{and}
\quad
h[n]=u[n-4]-u[n-15]
$$
By using distributivity property of the discrete-time convolution, we get :
\begin{align*}
(x*h)[n]&=(u[n-3]-u[n-8])*(u[n-4]-u[n-15])\\ \\
&=\underbrace{u[n-3]*u[n-4]}_{(y_{1}[n])} -\underbrace{u[n-3]*u[n-15]}_{(y_{2}[n])}-\underbrace{u[n-8]*u[n-4]}_{(y_{3}[n])}+\underbrace{u[n-8]*u[n-15]}_{(y_{4}[n])}
\end{align*}
Therefore we have four convolutions for which we can see that :
$y_{1}[n]$ = $\begin{cases}\displaystyle\sum_{3}^{n-4}1&\text{if $7\leq n \leq 11$}\\0&\text{if otherwise}\end{cases}=\begin{cases}n-6&\text{if $7\leq n \leq 11$}\\0&\text{if otherwise}\end{cases}$
$y_{2}[n] = \begin{cases}\displaystyle\sum_{3}^{n-15}1&\text{if $18\leq n \leq 33$}\\0&\text{if otherwise}\end{cases}=\begin{cases}n-17&\text{if $18\leq n \leq 33$}\\0&\text{if otherwise}\end{cases}$
$y_{3}[n] =\begin{cases}\displaystyle\sum_{8}^{n-4}1&\text{if $12\leq n \leq 16$}\\0&\text{if otherwise}\end{cases}=\begin{cases}n-11&\text{if $12\leq n \leq 16$}\\0&\text{if otherwise}\end{cases}$
$y_{4}[n]=\begin{cases}\displaystyle\sum_{8}^{n-15}1&\text{if $23\leq n \leq 38$}\\0&\text{if otherwise}\end{cases}=\begin{cases}n-22&\text{if $23\leq n \leq 38$}\\0&\text{if otherwise}\end{cases}$
I must note that I am afraid I don't know if the upper bound in the sum of $(1)$, $(2)$, $(3)$, and $(4)$ are right or wrong. I think what's left to do is compute :
$$
y[n]=y_{1}[n]-y_{2}[n]-y_{3}[n]+y_{4}[n]
$$
I hope someone can show me how from here we can prove that the answer is :
$$
y[n]=\begin{cases}
n-6&\text{if $7\leq n\leq 11$}\\
6&\text{if $12\leq n\leq 18$}\\
24-n&\text{if $19\leq n \leq 23$}\\
0&\text{if otherwise}
\end{cases}
$$


texts:
I have implemented a moving average, similar to the Hogenauer Filter, with a reduced number of computation operations. I expect the expected error to behave as the random walk and its STD to be of order $\sqrt{n} *\varepsilon$, where $n$ is the filtered vector length. Somehow I get one order smaller. What I am missing and is there a way to reduce the error?
len = 1000;
windLen = 11;
normCoeff = 1 / windLen;
q = nan(len, 1);

for a = 1:len
    x = randn(10^5, 1);

    xRef = movmean(x, windLen);
    xRef(1:(windLen - 1) / 2 + 1) = [];

    varState = 0;
    xTest = nan(size(x));
    x = [zeros(windLen, 1); x];
    for ind=1:length(x) - windLen
        varState = varState + x(windLen + ind) - x(ind);
        xTest(ind) = varState * normCoeff;
    end

    xTest(1:windLen) = [];
    xRef(length(xTest) + 1:end) = [];
    q(a) = xTest(end) - xRef(end);
end
disp(std(q))

edit
Following the suggestion by @Dan Boschen, I attached the comparison of this method vs Hogenauer Filter and for some reason, the attached method is faster (twice).
comment: Please mind that the external loop is just for the improvement of the comparison and not actually required for the computation.
clc
clear
windLen = 11;
testLen = 10^4;
normCoeff = 1 / windLen;
xBuff = zeros(windLen, 1);
x = randn(10^4, 1);

tic
for a = 1:testLen
    varState = 0;
    y = nan(size(x));
    xBuff(windLen + 1:windLen + length(x)) = x;
    for ind=1:length(x)
        varState = varState + xBuff(windLen + ind) - xBuff(ind);
        y(ind) = varState * normCoeff;
    end
end
toc

tic
for a = 1:testLen
    y2 = filter([1 0 0 0 0 0 0 0 0 0 0 -1], [11 -11], x);
end
toc
plot(y - y2)

The error accumulation appears here too.


texts:
I asked this on Math SE but received no replies. I hope this is a relevant forum to ask in.

I would like to analyze the spectrum of the following function:
$$f(t)=\cos(t\cdot a(1+b\cos(ct)))$$
with $a,b,c$ some real parameters. As the function is non periodic (perhaps unless $a,b,c$ are fine-tuned), the most natural approach is to Fourier transform it. But doing it by hand is quite horrible, Wolfram Alpha aso struggles, and using Mathematica I get
$$-ia\cdot\cos(\sqrt{2\pi})\cdot \delta'[\omega] - abc\cdot\cos^2(\sqrt{2\pi})\cdot \delta''[\omega]$$
I'm only interested in the real part of the result. I'm now trying to make sense of what I got.

1) Does the Mathematica result seem reasonable to you?

2) Looking at it from a generator function perspective (as in here), it seems there is a sharp peak centered at zero, accompanied by more minor peaks at $\pm\Omega$ for some $\Omega>0$. I find this result surprising, as I expected that the spectrum would be centered at $a$ as the "main" frequency followed by some smaller harmonics or some sort of $sinc$-like function.

3) Following a much simpler argument than the full Fourier transform, I try to think of it as a function that has a frequency that oscillates in "time" ($t$). Then maybe as an analogy I can ask myself which frequencies are most probable. In that case, I can think of the $\cos(ct)$ term as the cosine of a uniform random variable $ct:= Y\sim Uni(-\pi,\pi)$, and so the probability density function is actually U-shaped centered at $a$ and given by a function in the general form of $$f_Y(x)=\frac{1}{\pi\sqrt{1-x^2}}$$(normalized etc.)

But these two analyses are far off from one another, and I would like to understand what is going on.

A counter-argument to my probabilistic interpretation is that the term that determines the intensity of a frequency component is the coefficients of the Fourier decomposition (if it were periodic), and here this is not a linear combination of frequencies with different probabilities, but something else. Still, I need some help explaining the results. Thanks!


texts: Trying to understand the gold-rabiner pitch detection algorithm

texts:
I am currently in the process of taking on implementing the Gold-Rabiner pitch tracking algorithm, and I came across what seemed like a good explanation of it in David Meehan's "A Code Excited Linear Predictive Coder".

However, I am a bit confused with M5, which states:  "the magnitude of the difference in amplitude between the current valley amplitude and the previous peak amplitude"--  Where M3 is depicted as the "previous peak", yet the arrow is going up to where M1 is, which is the "current peak"-- so either the instruction is wrong or the diagram is wrong-- or I am misunderstanding this entirely.

Also, a side question is, because this is talking about magnitude, this means that the arrow pointing downward at M1, is pointing at 0, correct?


texts: Is it possible to replace an integrator system with an equivalent differentiator?

texts:
I have a system whose input-output relation is as follows

$$y(t)=x(t)+\int_{-\infty }^{t} x(\tau) \,\mathrm d \tau$$

Can I create an equivalent system by using differentiators rather than integrators?

I think something like taking derivative of both sides of the equation but improper integral makes it hard. Is there any nice way to convert this system? Thanks in advance.


texts:
If we let :
$$
x(t)=\begin{cases}
1&\text{if $0<t<1$}\\
0&\text{if otherwise}
\end{cases}
$$
and
$$
h(t)=x(t/a)=\begin{cases}
1&\text{if $0<t<a$}\\
0&\text{if otherwise}\end{cases}
$$
where $0<a\leq 1$
I wish to find convolution of $(x*h)(t)$ without graphing. We have that :
$$
\max(0,t-a)<\tau<\min(1,t)
$$

Is the following assertion correct?
$$
\max(0,t-a)=\begin{cases}
0&\text{if $0<t<a$}\\
t-a&\text{if $a<t<1$}
\end{cases}
\quad\text{and}\quad
\min(t,1)=\begin{cases}
1&\text{if $1<t<1+a$}\\
t&\text{if $0<t<1$}
\end{cases}
$$


texts: Amplitude value after a Discrete-Time integration calculation?

texts:
I am using Matlab block Discrete-Time
Integrator on Simulink and I would like to know in advance, what  output of that block regarding to its amplitude and the current sampling frequency would be.

For instance, I am setting as input  a $$ {\tt Reference} (t) = A \sin(\omega t)$$ where $ A = 0.4$ and $\omega = 2\pi\cdot 2KHz$, the outcome of  Discrete-Time
Integrator is my ${\tt Reference} (t)$ integrated a cosine function as expected  $$ve[n] = -B \cos[\Omega n]$$ with $B=2.55*10^4$

However, as shown in the picture  the amplitude $B$ is quite large, which puzzles me. I am aware that the integrator scales the output with $K T_s$, being my $f_s = 200MHz$. So I was wondering whether or not there is an expression that allows me know that scale factor in the amplitude that is happening after the integration, I am aware that is related to my sampling time but I cannot figure out how or what is.


texts: Unambiguously detecting (not correcting) between $t+1$ and $2t$ errors in a Reed-Solomon message

texts:
I know that Reed-Solomon codes of the form $(N,N-2t)$ can correct up to $t$ errors, and detect up to $2t$ errors, but I've only seen the correction procedure: Berlekamp-Massey or Euclidean to get the error-locating polynomial $\Lambda(x)$, then the Chien search to find its roots, then the Forney algorithm to find the error coefficients.

If I am told that there are at most $2t$ errors, is there a way to unambiguously tell how many errors $e$ there are? Or at least distinguish between the case of $e \le t$ (correctable errors) and $t < e \le 2t$? Is it just a matter of determining the degree of $\Lambda(x)$? Somehow that doesn't seem right, because then finding the roots of $\Lambda(x)$ would tell me where they were, and then I could correct them.

I'm confused, because both the Berlekamp-Massey and Euclidean algorithms appear to find polynomials of degree $2t$ or less -- for Berlekamp-Massey it's mentioned in Blahut's book, and for the Euclidean algorithm it's $\Lambda(x)S(x) + Q(x)x^{2t} = \Omega(x)$ for known syndrome polynomial $S(x)$.


texts:
I need to do a FFT on an image for noise reduction, but the problem is that I do not need the complete image, but only a circle in the middle. The borders are a fixed rig, thus I am not interested in what it displays, but it has an impact on the result of the FFT.

Is there any way to just cut out a circular part and use that for the FFT? Note that if I use black background, the edge between background, and image data will have quite an impact.


texts: Transient cancellation in non causal filtering

texts:
I am using the exponential moving average filter (IIR):
$$y(n)=\alpha{}\cdot{}x(n)+(1-\alpha{})y(n-1)$$
in a foward-backward filtering scheme, in order to eliminate completely any phase response phenomena. But I do not know how to eliminate the startup transients, ie. the transients from the initial state of the filter before and after the duration of the data.

Can you help me please?

I have in mind the Gustafsson's paper, but I found it a bit messy...


texts:
I have a transfer function coefficients as following:
[b,a] = ellip(4,.2,40,[.41 .47]);

I'd like to factor this system into its PFE using [r,p,k] = residue(b,a). Then, I'd like to implement a filter which uses a parallel combination of second-order subsections. To do this, I need to combine each complex pole with its complex conjugate so that the overall second order subsections will have real valued coefficients. Hence; I selected four pairs of first order terms and recombined them into second order subsections using "residue" function. How can I store the coefficients of the resulting second-order subsections in the matrices c and d so that each row corresponds to one second order system?

I tried to do this but I don't know whether it's correct or not?

[c1 d1]=residue([r(1) r(2)],[p(1) p(2)],0);
[c2 d2]=residue([r(3) r(4)],[p(3) p(4)],0);
[c3 d3]=residue([r(5) r(6)],[p(5) p(6)],0);
[c4 d4]=residue([r(7) r(8)],[p(7) p(8)],k);

[c d]=[c1 c2 c3 c4 d1 d2 d3 d4]

Can anybody please give me an idea? Any help would be appreciated


texts:
I'm using Digital Signal Processing Principles, Algorithm and Applications 4th Edition written by Proakis and Manolakis. In Chapter 3, subtopic 3.3.2 it mentions the term "causal signal" which I can't find any definition. I know what a causal system means: a system that output depend entirely on past and present input, not future input but here it's causal signal, not system

So what does causal signal definition?

I tried to look at the Index and it says that the term "Causal Signal" is mentioned in page 85 but in fact, I find nothing in page 85 mentioning it


texts:
Given a signal you are asked to perform baseband QPSK (not on a
  carrier).

Does this mean the following?

In baseband QPSK the (possible) constellation points are $$\sqrt {E_s} e^{jn \pi / 2} \qquad E_s = A^2 T_s$$

and the frequency response has delta functions placed in the origin

In QPSK modulation on carrier the (possible) constellation points are $$\sqrt {E_s} e^{jn \pi / 2} \qquad E_s = \frac {A^2 T_s} 2$$

and the frequency response has delta functions around the carrier frequency $\pm f_c$

Am I correct?


texts: How to prove $\int_{-\infty}^{\infty}x^2(t)dt=\int_{-\infty}^{\infty}\hat{x}^2(t)dt$ where $\hat{x}(t)$ is the Hilbert Transformation of $x(t)$

texts: $\hat{x}(t)$ is the Hilbert Transform of $x(t)$. Can anybody help me proving the above.


texts: Generate a nonuniform illumination (bias field) 2D image

texts:
I am creating a nonuniform illumination in MATLAB R2015. The nonuniform illumination (known as bias field) is defined as

Nonuniform Illumination (intensity inhomogeneity) which manifests itself as slow intensity variations in the same class over the image domain.
The non-linear degree of intensity inhomogeneity is indicated by the range of values of the bias field in the interval [1 − α, 1 + α] with α > 0

I used the below MATLAB code to generate the intensity inhomogeneity image in case of α = 0.2:
cols=256;rows=256;
alpha=0.2;
x = linspace(1-alpha,1+alpha,cols);
y = linspace(1-alpha,1+alpha,row2);
% Transform into a grid
[X,Y] = meshgrid(x,y);
bias=X;
imshow(bias,[])

However, I am confused about the term 'non-linear'. It means that the value in the range [1 − α, 1 + α] needs to be slowly changed in a non-linear manner. In the above code, it seems that it is linearly changing. Is this correct? Could you suggest me a method to create an image based on the definition?
This is an expected example which I need to obtain

or


texts: Cyclic prefix length in baseband OFDM with ZF equalizer

texts:
My question is about baseband OFDM. I try to plot the BER according to the SNR. I noticed that changing cyclic prefix length does not affect the performances. I.e., following the code below, gi = n_sp/1024 or n_sp/4 or even 0, gives the same performance (always the same BER curve). Is there a mistake somewhere?

NB: At the receiver I use ZF-equalizer. The full Matlab code is:

clear
close all
clc

n_subc = 2^11;                          % # of subcarriers
n_ofdm_sym = 10;                        % # of OFDM symbols
n_data_sym = n_subc*n_ofdm_sym;         % # of data symbols to transmit
gi = n_subc/4;                          % CP length

dqpskmod  =  comm.DQPSKModulator;
dqpskdemod = comm.DQPSKDemodulator;

h = [.7*exp(1i*pi/4) .4*exp(1i*pi/3) .2*exp(1i*pi/6) .1*exp(1i*pi/5)]; % 4-taps channel

SNR = -10:2:30;

for i = 1:length(SNR)
    for j = 1:100                                % Monte Carlo iterations

    d_symb = randi([0 3],1,n_data_sym);
    x_mod = step(dqpskmod, d_symb').';           % DQPSK modulated symbols

    x_mod_s2p = reshape(x_mod,n_subc,[]);        % serial to parallel

    x_ifft = ifft(x_mod_s2p,n_subc);

    x_cp = [x_ifft(end-gi+1:end,:); x_ifft];     % add CP

    x_ifft_p2s = reshape(x_cp,1,[]);             % parallel to serial

    % --------------------> Channel mixing and adding noise

    y = conv(x_ifft_p2s,h,'same');                          % Linear convolution
    % y = ifft(fft(h,length(x_ifft_p2s)).*fft(x_ifft_p2s)); % OR circular convolution
    y = awgn(y,SNR(i),'measured');                          % add noise

    % --------------------> Receiver

    y_s2p = reshape(y,n_subc+gi,[]);             % S2P

    y_cp = y_s2p(gi+1:end,:);                    % Removing CP

    y_fft = fft(y_cp,n_subc);                    % FFT

    y_p2s = reshape(y_fft,1,[]);                 % P2S

    pilotX = x_mod(1:4:end);                     % Pilot extraction from the DQPSK transmitted data
    pilotY = y_p2s(1:4:end);                     % Pilot extraction from the DQPSK received data
    H_est = pilotY./pilotX;                      % Channel estimation
    H_est_interp = interp1(1:4:n_data_sym,H_est,1:n_data_sym); % Linear interpolation
    H_est_interp(end) = H_est_interp(end-1);     % To avoid the last NaN

    X_est = y_p2s./H_est_interp;                 % ZF equalizer

    y_demod = step(dqpskdemod,X_est.')';         % DQPSK demapping

    d_symb(1:4:end)  = [];                       % Pilot suppression from the original data
    y_demod(1:4:end) = [];                       % Pilot suppression from the received data

    [~,ber(j)] = biterr(d_symb,y_demod);

    end
    BER(i)=mean(ber)
end

figure
semilogy(SNR,BER)

The figure below gives the results for gi = n_sp/1024, n_sp/4, and 0:

Additional question: why BER saturates for SNR above 22dB? There is no saturation for QAM modulation, is it specific to the DQPSK?


texts: Are EEG signals periodic or a-periodic, and if they are a-periodic, how can we assign frequencies to them.


texts: What is the relationship between angular frequency and normalized angular frequency

texts:
This is a slide from my lecture notes:

My professor used the following words " We denote digital frequencies with capital letters and analogue frequencies with lower case letters"

The problem I have is $$\Omega=\omega T_s$$

Wherever I am looking it is actually the other way round i.e.

$$\omega = \Omega T_S$$

For example  here

Is my proffesor wrong?


texts: How can we prove that Spectral density of In-Phase component and Quad-Phase component of a sample function of narrow-band noise have same definition?

texts:
Suppose we have a narrow-band noise $N(t)$. We analyse this to In-Phase and Quad phase components by the equation:

$$N(t) = N_I(t)  \cos(2\pi f_c t) - N_Q(t) \sin(2\pi f_c t)$$

Now a book (by Simon Haykin) states that:

$$ S_{N_I}(f) = S_{N_Q}(f) =
\begin{cases}
 S_N(f-f_c) + S_N(f+f_c) & -B\le f \le B \\
 0 & \text{otherwise}
\end{cases}
$$

How do we prove this?

This is the Image of Power Spectral Density of Narrowband Noise.


texts: I am reading about CFO estimation in OFDM system, assume for example that we estimate it based on CP as mentioned HERE and the estimated CFO for the received signal $r(t)$ is $β$. My question, how now can we compensate that CFO on the time-domain signal $r(t)$ ?


texts: Extracting direction of dominant sinusoids from an image in terms of angles

texts:
Here is a sample image.

What I am looking to do is compute the dominant "wrinkle" directions.
The dominant wrinkle directions can be seen in the following FFT, roughly, by drawing a line through the highest valued set of points in the spectrum.

I have been told that a potential way to find this direction would be to look at the auto-correlation scipy.signal.correlate2d.
I have not worked with this tool before, so I want to make sure my idea is correct.
The way that I think this will work is by using the FFT data to find the largest components of the signal.
Say for example, I look at the largest two components.
I can construct two images, one for each of the largest components, lets call the images dominant1 and dominant2.

For each image, I could then look at:

scipy.signal.correlate2d(image,dominant1)
scipy.signal.correlate2d(image,dominant2)

I believe that this will return a matrix that shows where the features are in the image that most agree with the signals dominant1 and dominant2.

If this interpretation is correct, I am still unsure of how to extract the dominant_angle1 and dominant_angle2 say. Suppose that the horizontal direction in the image as shown above is x and the vertical is y. I would like dominant_angle1/2 to be the angle off of the x-axis.

Any comments, suggestions for better approaches, will be appreciated.


texts: Conceptual question on FFT/IFFT (IFFT existence)

texts:
I was reading  "Discrete and continuous Fourier transforms: analysis, applications and fast algorithms" written by E.Chu and, at some point, I found something that I could not completly understand.

  If $x(t)$ is not continuous, then it is possible that 2 different functions share the same transform $X(f)$. Therefore, the function $x(t)$ cannot be uniquely determined by inverse trasforming $X(f)$ unless it is continuous or it satisfies 2 more conditions:

  1- $x(t)$ has only a finite number of maxima and minima on any finite interval

  2- $x(t)$ has on any finite interval at most a finite number of discontinuities, each of which is a jump discontinuity.

  In the latter case, the inverse Fourier transform of $X(f)$ produces $\hat{x}(t)$, which agrees with $x(t)$ at every $t$ at which $x(t)$ is continuous, and $\hat{x}(t)$ equals the average of the left-hand limit $x(t^-_{\alpha})$ and the right-hand limit $x(t^+_{\alpha})$ at every $t_{\alpha}$ at which a jump discontinuity occurs.

So, here my doubts.

First of all, I am not sure I got the meaning of the statements related to the 2 conditions listed. This is probably due to the fact that I am familiar only with signals coming from structural tests (mainly accelerations), therefore I think it could be really helpful if someone of you can provide me some example I can focus on to understand the concept (for example, how can be possible to have infinite maxima and minima in a signal? And what are the discontinuities it is talking about?).

Second, measured signals (for example, the aformentioned acceleration measurements coming from a test) are not continuous (since they are sampled at a given sampling frequency), so if I want to apply the IFFT algorithm to retrieve the time signal does it mean I have always to check if the 2 other conditions are satisfied? And how can I check them? Visual inspection? Or do I have to look for some specific features the signal must have?

Are there classes of signals which always satisfy the 2 conditions?


texts:
can someone define/explain briefly (or with some detail-> up to you) the following terms:

Block Processing? and what are the differences between it and Sample Processing?
(removed)
What is a 3D Audio System and what is relation between it & Head Related Transfer Function (HRTF)? and how can HRTF be applied in Headphones?
(removed)
The relationship between bit rate, sample rate and bit depth? Is there a formula that combines all terms?


texts: How to prepare different input size in CNN

texts:
I'm doing a research about personality identification based on their signature using CNN method, however the learning feature for the personality traits have a different input size. I understand that in CNN we have to make the same size for the input, but can I make the input size in different size?

Here is the journal that I use for my references for the learning feature.